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Preface 


Like  most  of  my  peers,  I  have  long  been  interested  in  space  and  most  anything  related  to  space 
exploration  and  utilization.  After  being  assigned  to  first  Air  Force  Space  Command  and  then  United 
States  Space  Command,  I  found  that  I  enjoyed  studying  how  objects  travel  through  spa/x  and  eq)ecially 
how  to  track  and  predict  their  movements  around  the  earth.  So  it  was  very  satisfying  to  be  able  to  work  on 
a  project  that  has  both  theoretical  and  practical  applications. 

As  an  orbital  analyst  at  the  Space  Surveillance  Center  in  Cheyenne  Mountain,  1  routinely  tradced 
objects  and  updated  their  predicted  positions  but  I  could  not  have  explained  how  the  mathematics  woiiced. 
Through  the  classes  in  Orbital  Dynamics  taught  here  at  AFTT,  1  learned  the  theory  of  how  orbits  could  be 
estimated  using  observations.  But  it  was  the  research  for  this  thesis,  the  development  of  programs  to 
simulate  the  orbit  of  a  satellite  and  to  generate  radar  observations  and  the  use  of  these  observations  to 
estimate  the  state  of  the  satellite,  that  fully  developed  the  level  of  understanding  that  I  desired.  This 
project  has  been  both  difficult  and  challenging,  but  few  worthwhile  endeavors  are  not. 

1  could  not  have  completed  this  woilc  without  the  help  and  encouragement  of  my  teachers,  friends 
and  ftunily.  I  would  like  to  thank  Dr.  Wiesel  for  his  patience  and  understanding  as  well  as  for  his  insight 
and  direction  when  things  got  tough.  1  would  also  like  to  thank  my  fellow  classmates,  eq)ecially  Ctqrtain 
Tony  Nash,  for  answering  my  many  questions  and  for  simply  listening  to  my  ideas.  And  last,  but 
certainly  not  least,  my  wife  who  suffered  through  the  entire  process  with  me. 
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Abstract 


This  study  investigated  the  use  of  a  least  squares  differential  corrector  to  compress  the 
information  contained  in  a  single  radar  track  of  a  Near  Earth  satellite  into  a  six  element.  Earth  Centered 
Inertial  state  vector  and  associated  covariance  matrix.  Observations  were  generated  using  a  truth  model 
program  based  on  two-body,  J2  geopotential,  and  atmospheric  drag  dynamics  and  consisted  of  simulated 
range,  azimuth  and  elevation.  Random  Gaussian  noise  was  added  to  the  data  to  simulate  the  randmn 
errors  seen  in  real  world  observations.  The  orbital  dynamics  used  in  the  estimator  consisted  of  two-body 
and  J2  geq)otential  effects  which  were  approximated  using  a  Taylor  Series  expansion  in  time.  Analysis 
was  conducted  using  two  circular,  inclined  low  earth  orbits.  Estimator  performance  was  evaluated  as  a 
function  of  maximum  pass  elevation,  track  duration,  data  rate  and  the  highest  order  of  approximation  of 
the  equations  of  state.  The  estimator  successfully  extracted  all  of  the  pertinent  information  from  the  data 
with  accuracy  being  dependent  on  the  number  of  data  points  available.  A  correlation  between  the  requited 
order  of  approximation  and  maximum  expected  track  length  was  also  noted. 
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COMPRESSION  OF  A  RADaR  TRACK  OF  A  NEAR  EARTH  SATELLITE  INTO  AN 
EARTH  CENTERED  INERTIAL  STATE  VECTOR  USING  LEAST  SQUARES 

DIFFERENTIAL  CORRECTION 


I.  Background  and  Statement  of  the  Problem 


1. 1  Introduction 

The  Space  Surveillance  Center  (SSC),  located  in  Chejvnne  Mountain  Air  Force  Base.  Colorado,  is 
the  world's  premier  space  tracking  organization.  Through  the  use  of  a  world  wide  network  of  ground 
based  radars  and  electro-optical  cameras  known  as  the  Space  SurveillaiKe  Network  (SSN),  the  SSC  trades 
and  maintains  orbital  predictions  for  all  man-made  objects  in  orbit  about  the  earth  with  a  radar  cross 
section  of  at  least  ten  square  centimeters.  Satellites  with  a  period  less  than  225  minutes  are  defined  as 
Near  Earth  satellites.  Almost  all  of  the  data  used  to  trade  near  earth  satellites  comes  from  ground  based 
qrace  surveillance  and  ballistic  missile  early  warning  (BMEWS)  radars. 

When  these  radars  track  a  satellite  th^  produce  data  called  observations  which  are  sent  to  the  SSC 
over  a  teleconununication  link.  In  their  most  basic  form  this  data  consists  of  the  time,  elevatum,  azimuth 
and  range  of  the  satellite  as  observed  from  the  site.  A  typical  track  for  a  near  earth  ol:^  might  span 
from  five  to  twenty-five  minutes.  Since  some  of  these  radars  can  rtteasure  the  position  of  a  satellite  as  first 
as  100  tirr^  each  second,  this  could  result  in  30,000  to  150,000  observations  per  trade  (7:  156). 
Obviously,  transmitting  this  much  data  to  the  SSC  every  time  a  satellite  is  tracked  is  not  feasible.  Instead, 
the  radar  sends  only  a  representative  sampling  of  this  data  to  the  SSC.  Typically,  this  consists  of  three 
observations  positioned  at  the  beginning,  middle  and  end  of  the  track.  For  very  low  priority  satellites 
perhaps  only  one  observation  will  be  sent  to  the  SSC. 
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The  SSC  combines  this  new  data  along  with  observations  from  other  sensors  using  batch  and/or 
sequoitial  processing  to  update  the  satellite's  orbital  elements  or  elset.  Every  catalogued  satellite  still  in 
orbit  has  its  own  elset  which  is  continuously  maintained  using  the  observations  sent  in  by  the  SSN.  These 
elsets  are  then  used  to  predict  satellite  positions  in  the  near  future  using  both  general  and  special 
perturbation  methods.  One  important  use  for  these  updated  elsets  is  to  allow  sensors  to  acquire  and  track 
satellites.  Maintaining  an  accurate  and  timely  database  is  also  important  in  avoiding  satellite  collisions 
and  providing  warning  of  over  flights  of  foreign  intelligence  satellites. 

With  over  7,000  objects  in  orbit  and  more  than  twenty  six  sensor  sites  in  the  SSN,  the  SSC  typically 
receives  30,000  to  80,000  observations  per  day.  This  is  a  large  aiiKMmt  of  data  to  be  transmitted  to  aitd 
processed  by  the  SSC  computer.  As  space  continues  to  increase  in  importance  for  both  military  and 
conunercial  purposes,  the  number  of  objects  in  earth  orbit  will  iiKiease.  This  will  result  in  additional 
woric  loads  for  both  the  sensors  making  up  the  SSN,  for  the  SSC  computer  and  persoiutel,  artd  for  the 
communications  links  that  tie  them  together. 

If  a  way  could  be  found  to  more  efificiently  use  the  information  gathered  by  the  sensors  in  the  SSN  it 
might  be  possible  to  alleviate  some  of  the  heavy  demands  that  can  be  expected.  This  would  allow  existing 
resources  to  handle  the  job  without  expensive  investments  to  upgrade  facilities  and  equipment.  One  way 
to  do  this  would  be  to  compress  all  of  the  data  into  a  single  state  vector  and  associated  covariance  matrix. 
This  would  allow  all  of  the  information  contained  in  the  track  to  be  sent  to  the  SSC  instead  of  a  very  small 
sartq>le.  Not  only  would  this  provide  the  SSC  with  better  data  which  could  reduce  the  frequency  of 
required  updates  or  sensor  tasking,  but  it  distributes  some  of  the  processing  load  out  to  the  sensor  sites. 

1.2  Problem  Statement 

In  order  to  increase  the  capacity  of  the  SSC  and  SSN  to  track  and  maintain  elsets  on  the  growing 
poinilation  of  near  earth  space  objects,  it  is  necessary  to  develop  a  method  of  compressing  all  of  the 
information  contained  in  a  sensor  track  into  an  easily  transmitted  and  utilized  form.  The  m^hod  must  be 
easily  implemented  and  not  require  significant  equipment  upgrades.  A  possible  solution  is  to  have  each 
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smsor  site  convert  all  of  its  observations  into  a  state  vector  and  associated  covariance  matrix  using  a  least 
squares  batch  differential  correction  based  on  two-body  dynamics  and  possibly  J2  effects. 

1.3  Research  Objectives 

The  first  objective  is  to  develop  a  truth  model  to  generate  simulated  ground  sensor  site  observations. 
The  truth  model  should  model  satellite  motion  using  two-body  orbital  dynamics,  the  J2  zonal  harmonic, 
and  air  drag  since  these  are  the  three  dominant  factors  in  the  orbital  motion  of  near  earth  satellites. 

Output  from  the  truth  model  will  consist  of  accurate  time,  range,  azimuth,  and  elevation  values  to  which 
Gaussian  random  noise  will  be  added  to  simulate  real  wortd  observations. 

The  second  objective  is  to  create  a  differential  corrector  based  on  least  squares  estimation  (qrerating  in 
a  batch  mode.  The  dynamics  model  used  in  the  estimator  will  be  based  on  a  Taylor  Series  expansion  of 
the  equations  of  motion  for  an  object  in  orbit  about  the  earth.  The  state  transition  matrix  used  in  the 
estimator  will  be  derived  from  these  same  Taylor  Series  equations. 

The  third  objective  is  to  investigate  the  effectiveness  of  the  differential  corrector  in  reducing  the 
simulated  observations  from  the  truth  model  into  a  single  state  vector  and  associated  covariance  matrix. 
Specifically,  it  will  be  necessary  to  determine  how  many  terms  are  needed  in  the  Taylor  Series 
approxirruition  of  the  equations  of  motion.  Also  it  must  be  determined  if  the  equations  of  motion  must 
include  the  J2  zonal  harmonics  or  if  the  two-body  dynamics  alone  are  sufficient.  And  finally,  the 
correlation  between  pass  time  (or  maximum  elevation)  and  the  effectiveness  of  the  estimator  will  be 
studied. 

1.4  Research  Questions 

Least  squares  estimation  is  a  commonly  used  method  for  performing  differential  correction  on  the 
orbital  elements  or  state  of  satellites.  It  is  especially  effective  when  used  in  a  batch  mode  where  sufficient 
observations  exist  to  span  a  significant  portion  of  the  object's  orbit.  In  the  case  of  a  single  track  by  a 
ground  based  sensor  site  of  a  near  earth  satellite,  only  a  srruill  part  of  the  orbit  will  be  observed.  Two 
factors  that  will  influence  the  number  of  observations  will  be  the  objects  maximum  elevation  relative  to 
the  site  and  its  altitude. 
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This  leads  to  the  following  questions: 

1.  Does  it  matter  how  good  the  dynamics  model  used  by  the  estimator  is? 

a.  Are  the  effects  of  the  J2  zonal  harmonic  or  of  atmo^heric  drag  large  enough  to  be  significant? 

b.  How  many  terms  must  be  included  in  the  Taylor  Series  expansion  of  the  equations  of  state? 

2.  What  is  the  minimum  number  of  observations  required  by  the  estimator  to  effectively  conqrtess  the 
data  into  a  single  state  veaor  and  covariance  matrix? 

3.  Can  this  be  related  to  the  maximum  elevation  of  the  pass  and/or  the  satellites  altitude  during  the 
pass?  If  so,  what  is  this  relationship? 

I.S  Assumptions 

Although  there  exist  numerous  perturbation  sources,  only  those  perturbations  caused  I9  the  J2  zonal 
effect  and  atmospheric  drag  are  considered  in  the  truth  model.  This  is  based  on  the  reasonable 
assumption  that  the  other  perturbations  are  too  small  to  be  noticeable  -  especially  during  the  short  track 
times  for  near  earth  objects. 

Normally,  several  d^s  worth  of  data  are  required  to  estimate  the  effect  of  drag  on  an  orbit.  The  SSC 
typically  uses  observations  spaced  over  a  period  of  3-7  days  to  determine  the  drag  experienced  by  a  Near 
Earth  satellite.  Because  this  study  lodes  at  reducing  the  data  from  a  single  track,  insufficient  data  is 
available  for  estimating  the  effect  of  drag.  Therefore  the  estimator  will  use  only  two-body  and  J2  effects. 

The  truth  model  generates  accurate  ephemeris  information  for  the  input  satellite  reference  trajectory. 
The  data  is  assumed  to  be  perfect  and  is  made  more  realistic  by  applying  Gaussian  random  noise  to  the 
data. 

There  are  errors  associated  with  the  dynamics,  the  state  vector,  and  the  observations.  TIw  estimation 
routines  developed  during  this  study  assurtK  that  the  use  of  a  Taylor  Series  expansion  will  result  in  a 
(fynamics  model  that  is  close  enough  to  reality  and  that  the  major  source  of  uncertainty  comes  from  the 
observations. 
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1.6  Scope 

The  pr(qx>sed  study  models  the  effectiveness  in  a|q>lying  a  two-body  least  squaies  batch  estimator  to  a 
short  arc  of  observations  to  create  a  single  state  vector  and  associated  covariance  matrix.  The  main  efforts 
in  this  tl^is  are:  (1)  developing  a  truth  model  to  generate  simulated  observations  from  ground  based 
sensor  sites,  (2)  building  a  least  squares  differential  corrector  to  generate  an  earth  centered  inertial  (ECI) 
state  veaor  and  covariance  matrix,  and  (3)  analyzing  the  effectiveness  of  this  estimator  in  reducing  or 
compressing  the  data  for  near  earth  satellite  passes. 

1.7  Summary 

As  the  number  of  objects  in  earth  orbit  continues  to  increase  along  with  the  growing  irtqxHlance  of 
space  in  general,  the  demands  made  on  the  SSC  and  SSN  to  track  and  predict  the  motion  of  these  objects 
will  also  increase.  Instead  of  trending  large  sums  of  money  to  upgrade  facilities  and  equipment  it  may  be 
possible  to  multiply  the  effectiveness  of  current  capabilities  of  the  SSC  and  SSN. 

One  possible  solution  is  to  compress  all  of  the  information  from  a  single  sensor  track  so  that  it  can  be 
easily  transmitted  to  and  used  by  the  SSC.  This  could  be  done  I9  using  a  least  squares  batch  differential 
corrector  to  reduce  the  data  obtained  at  the  sensor  site  into  a  single  state  vector  and  associated  confidence 
(covariance)  matrix.  This  would  essentially  extract  more  usable  data  from  each  sensor  trade.  It  is  hoped 
that  this  would  allow  for  more  accurate  elsets  to  be  maintained  at  the  SSC.  This  in  turn  could  reduce  the 
frequency  with  which  elsets  must  be  updated  and  would  essentially  preprocess  the  observations,  thus 
redudng  the  computational  load  on  the  SSC  computer.  Additionally,  it  might  be  possible  to  reduce  the 
number  of  tracks  or  the  fraction  of  a  track  required  on  each  satellite,  therd^  allowing  the  SSN  to  trade 
more  objects  with  existing  or  reduced  resources. 


11.  Theoreticai  Background 


2.1  Introduction 

For  almost  as  long  as  man  has  been  on  earth,  he  has  looked  tq>  into  the  night  sky  and  wondered  about 
the  things  he  saw  there.  As  far  back  as  3000  B.C.,  Babylonian  priest-astrologers  were  observing  and 
recording  the  motions  of  the  moon,  sun  and  the  five  visible  planets.  As  time  passed,  ancient  man 
attempted  to  devise  explanations  for  how  the  planets  moved  but  it  was  not  until  the  early  1600's  that 
Kqrler  finally  developed  his  three  laws  which  correctly  described  the  motion  of  these  heavenly  bodies. 
However,  Kepler's  explanation  merely  stated  how  the  planets  (or  any  other  object  in  qmoe  for  that  matter) 
moved  -  it  did  not  answer  the  more  fundamental  question  of  why.  Although  (jalileo  had  a  few  ideas,  it 
was  Isaac  Newton  who  finally  answered  the  question  when  he  published  his  three  laws  of  motion  and  the 
law  of  universal  gravitation  (l;7S7-738). 

The  first  method  for  finding  the  orbit  of  an  object  from  observations  was  published  Newton  in  the 
Principia.  This  method  required  three  observations  (angles  only  since  radars  did  not  exist  in  1703)  and 
consisted  of  a  graphical  approach  using  successive  approximations.  The  first  analytical  method  for 
solving  this  same  problem  was  presented  by  Euler  in  his  Theory  of  the  Motion  of  Planets  and  Comets  in 
1744.  Euler  also  developed  a  relationship  which  related  ar^r  two  position  vectors  with  the  vector 
connecting  them  and  the  time  interval  between  the  two  vectors.  Further  refinements  to  Euler's  work  were 
produced  by  Lambert  and  Lagrange  but  it  was  not  until  Laplace  develqred  a  new  method  in  1780  that  it 
was  possible  to  process  more  than  three  observations  at  a  time  (2:31-33,  117-122).  By  1802,  Karl 
Frederick  Gauss  had  developed  yet  another  method  for  using  three  observations  to  d^rmine  the  orbital 
elements  of  an  object,  the  basic  elements  of  which  are  ^11  used  today. 

As  important  as  all  these  developments  were,  an  even  more  important  contribution  to  the  science  of 
orbit  determination  was  Gauss'  method  of  least  squares.  By  using  the  method  of  least  squares  it  was 
possible  to  use  more  than  three  observations  to  determine  the  orbit  of  an  object  to  fer  greater  precision 
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than  before  by  producing  orbital  elements  that  agreed  with  all  of  the  observations  used  rather  than  with 
just  two  or  three.  As  Gauss  himself  wrote: 

Finally,  as  all  our  observations,  on  account  of  the  imperfection  of  the  instruments  and  trf  the 
senses,  are  only  ai^roximations  to  the  truth,  an  orbit  based  only  on  the  six  absolutdy  necessary 
data  may  still  be  li^le  to  considerable  errors.  In  order  to  diminish  these  as  much  as  possiUe, 
and  thus  to  reach  the  greatest  precision  attainable,  no  other  method  will  be  given  except  to 
accumulate  the  greatest  number  of  the  most  perfect  observations,  and  to  adjust  the  elements,  not 
so  as  to  satisfy  this  or  that  set  of  observations  with  absolute  exactness,  but  so  as  to  agree  with  all 
in  the  best  possible  manner.  (3) 

2.2  The  Method  of  Least  Squares 

The  fundamental  elements  of  least  squares  are  built  on  four  basic  assumptions:  that  the  dynamics 
contain  no  error,  that  the  observations  do  contain  errors,  that  the  observational  errors  obey  a  Gaussian 
distribution,  and  that  the  Principle  of  Maximum  Likelihood  ai^lies.  The  last  item  comes  from  Gauss' 
groundbreaking  work  in  the  field  of  probability  theory  and  states  that  the  value  of  an  estimate  which 
maximizes  the  probability  of  obtaining  a  particular  set  of  data  is  the  true  value  (4:  22).  The  third 
assumption  is  supported  by  the  Central  Limit  Theorem,  also  develq)ed  by  Gauss,  which  states  that  the 
errors  in  the  observation  data  will  follow  a  normal  distribution  (see  Figure  1)  as  long  as  the  total  error  is 
the  product  of  many  small  error  sources  (4:  1 1). 


z 

Figure  1.  Gaussian  Distribution 
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The  assumption  that  the  errors  associated  with  the  observation  data  follow  a  normal  or  Gaussian 
distribution  is  a  key  element  in  the  development  of  estimation  theory.  Therefore  the  probability  of  an 
error  is  governed  by  the  Gaussian  probability  density  Junction: 


P(e)de  = 


1 


(1) 


where  P(e)  is  the  probability  that  an  error  e  will  occur.  This  implies  that  the  probability  of  an  error  is 
defined  by  the  area  under  the  curve  containing  the  error  in  question.  Thus  it  is  necessary  to  define  a  range 
of  errors,  swrh  as  e  and  e+Se,in  order  to  calculate  the  probability.  Eq  ( 1)  and  Figure  1  demonstrate 
several  important  features  of  a  Gaussian  distribution;  (1)  The  distribution  is  symmetric  about  the  true 
value  of  the  measurement.  (2)  The  width  of  the  curve  is  determined  by  the  standard  deviation  a.  (3)  The 
tmal  area  under  the  curve  is  one  because  the  function  is  normalized.  This  means  that  there  will  alw^  be 
some  error  in  a  measurement,  although  small  errors  are  more  likely  than  large  ones.  A  large  value  of  o 
results  in  a  broad  curve,  signifying  low  precision,  while  a  small  value  of  o  results  in  a  narrow  curve, 
signifying  high  precision. 

The  next  step  in  the  development  of  estimation  theory  irtvolves  the  estimation  qrerator; 

E(-)  =  (2) 

where  the  qien  slot  (-)  contains  the  quantity  to  be  estimated  and  f(x)  is  a  probability  density  function.  By 
applying  the  estimation  operator  three  special  results  are  obtained:  (1)  The  erqrected  value  of  a 
measurement  is  the  true  value:  ECx)  =  Xg  .  (2)  The  average  error  in  the  measurement  is  zero;  E(e  =x- 
Xg)  =  0  .  (3)  The  average  squared  error  is  a^;  E(e^)  =o^  (4:19-20). 

For  the  case  where  only  a  single  value  is  to  be  estimated,  such  as  the  length  of  a  rod,  the  probability  of 
obtaining  a  particular  set  of  independent  measurements  is  given  l^: 


P{x^,x^ 


,...Xfj  )  —  (27C) 


-NI2 


(3) 
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where  Oj  is  the  standard  deviation  associated  with  each  measurement.  A'  is  the  total  number  rtf' 
measurements,  Xj  is  the  value  of  each  measurement,  and  is  the  true  value.  Since  the  true  value  can 

never  be  known,  it  is  necessary  to  apply  the  Principle  of  Maximum  Likelihood  which  allows  the 
replacement  of  with  its  estimate  x  (4:22).  In  order  to  maximize  the  probability  of  the  data  with 
respect  to  the  estimate,  the  argument  of  the  exponential  must  be  minimized; 


=  0 


(4) 


giving 


I 


i=l 


(x,-X) 

^2 


=  0 


(5) 


It  is  from  Eq  (4)  that  this  estimation  method  derives  its  name  of  least  squares.  The  quantity  (x^-x)  is 
usually  called  the  residual,  r/ ,  which  can  be  thought  of  as  an  estimate  of  the  true  error.  Solving  for  the 
estimate  x  gives; 


2.S  Multi-Dimensional  Least  Squares 

The  preceding  section  showed  the  development  of  least  squares  estimation  theory  for  a  one 
dimensional  or  single  variable  case.  However,  most  problems  (such  as  orbit  estimation)  are  multi¬ 
dimensional  in  nature.  Fortunately,  the  results  from  section  2.3  can  be  extended  to  siqrport  multi¬ 
dimensional  problems.  To  start,  define  the  following  notation: 

=  iX^,Xj,X^,...,X^)  (7) 

and 
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0,1  0,2  ...  Oiff 

®21  ®22  •••  ®11 

_2  2  2 
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where  is  the  slate  or  variables  in  the  problem  and  P  is  the  covariance  matrix.  The  covarianoe  matrix 
contains  information  on  the  variance  associated  with  each  element  in  the  state  vector.  The  diagonal 
entries  are  called  variances  while  the  off-diagonal  entries  are  referred  to  as  covariances.  Using  this 
notation  the  Gaussian  probability  density  function  can  be  written  as; 

/(i)  =  |Pr'"  expf-i(i  - 1. )'  P-'  (I  - 1,, )  1  (9) 

where  the  variable  x  has  been  replaced  by  the  estimated  state  x ,  and  the  standard  deviation  squared, 
has  been  replaced  by  the  covarianoe  matrix  P.  As  before,  the  argument  of  the  exponential  is  the  quantity 
to  be  minimized,  yielding; 

(x-Xo)’'p-'(i-Xo)>0  (10) 

which  means  the  covariance  matrix  P  must  be  a  positive  semi-definite  matrix  (4;27). 

2.4  Linearizing  the  Problem 

The  problem  of  orbit  determination  is  fundamentally  a  nonlinear  one.  Not  only  are  the  dynamics 
highly  nonlinear  but  so  is  the  relationship  between  observations  and  the  satellite's  state.  However,  in 
order  to  make  any  progress  in  estimating  the  state  of  a  satellite  it  is  necessary  to  linearize  various  aq)ects 
of  the  problem. 

To  start,  the  dynamics  must  be  specified,  either  as  a  set  of  differential  equations 

^  =  (11) 

or  as  an  explicit  solution 
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x(/)  =  h(i(lo).0 


(12) 


which  gives  the  state  as  a  function  of  time  and  the  initial  condititms.  Since  the  dynamics  are 
detenninistic,  these  equations  can  be  linearized  as: 

5x(/)  =  <I>(l,/o)5*(/o)  (13) 

where  is  the  state  transition  matrix  relating  a  change  in  the  state  at  time  /  to  a  change  in  the  initial 
conditions  at  time  tp.  If  the  (lynamics  are  given  e)q)licitly  as  a  function  of  time  then  the  state  transition 
matrix  is  given  by: 

<*>(/,/o)  =  V.(^)h(x(/o),/)  (14) 

As  was  pointed  out  earlier,  the  true  state  Xg  will  always  be  unknown  and  since  the  whole  purpose  of  least 
squares  is  to  produce  an  estimate  x,  a  reyference  5/0/e  must  be  specified  about  which  the  dynamics 
can  be  linearized.  The  relationship  between  the  reference  state  and  the  estimate  is: 

x=x„^+8x  (15) 

where  5x  is  the  correction  needed  to  turn  the  reference  state  into  the  estimate.  Hopefully  the  initial  guess 
of  the  state  will  be  close  enough  to  the  true  state  to  allow  the  linearized  dynamics  to  be  valid.  By  making 
successive  corrections  to  the  reference  state,  it  should  be  possible  to  converge  to  the  correct  estimate  of  the 
true  state. 

The  relationship  between  the  observations  and  the  ^em  state  can  be  written  as 

z,(0  =  G(x(0./,)+e,  (16) 

where  Zj  is  the  observation  vector,  such  as  {range,  azimuth,  elevation},  G  is  the  function  which  converts 
the  state  variables  at  time  tj  to  the  observation  variables,  and  Cj  is  the  true  error  associated  with  the 
observation.  Since  this  relationship  is  also  usually  nonlimar  it  will  be  necessary  to  linearize  it  as  well: 
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(17) 


e  =  *-*o 

=  G(x,/)-G(»„,/) 

=  G(Xo  +5x,0-G(Io,0 


dx 


5x(0 


where  z  is  the  actual  observation  vector,  Zg  is  the  observation  vector  using  perfect  data,  x  is  the  observed 
state,  Zq  is  the  true  state,  and  5x  is  the  difference  between  the  true  and  observed  states.  Replacing  x  in  Eqs 
(16)  and  (17)  with  x^^f  allows  the  definition; 

«,).<,)  (1«) 


which  relates  the  error  in  the  state  to  the  error  in  the  reference  trajectory  (4:69). 

2.S  Correcting  the  Reference  State 

Just  as  before,  the  concept  of  finding  the  residual  instead  of  the  true  error  must  be  used 

r,  =2,-G(x,^(f,),/,)  (19) 


where  the  residual  vector  r  is  formed  from  the  observed  location  z  minus  the  predicted  location  G(x).  If 
the  reference  state  is  close  enough  to  the  true  state,  then  the  residual  vector  will  be  linearly  related  to 
changes  in  the  reference  state  (4:70); 

5r,  *H,6i(/,)=H,<I>(/,fo)6x(/o) 

=  TMto) 

and 

r/  «  +  6r,  =  )  (21) 

where  r*  is  the  residual  vector  after  the  reference  state  is  corrected  to  be  the  estimate,  and  rj  is  the  residual 
vector  associated  with  the  reference  state.  The  matrix  T  is  called  the  observation  matrix  artd  is  defined  as 
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T,  sH,<D(/,,/o) 


(22) 


Since  least  squafes  is  normally  used  with  multiple  observations  it  will  be  conveniem  for  notation 
purposes  to  define  three  items.  The  first  is  the  total  residual  vector 

r^^[r[  r[  ...  r^]  (23) 

where  each  rj  is  itself  a  vector  of  residuals.  The  second  is  the  observation  matrix 


■q,  0  ...  0 

0  0  ... 

where  each  Qi  is  itself  a  matrix  of  the  instrumental  covariance  associated  with  each  data  vector  Zj.  In  the 
form  shown,  each  data  vector  is  statistically  independent  and  therefore  the  off-diagonal  entries  are  zero. 

In  accordance  with  the  Principle  of  Maximum  Likelihood,  the  estimate  x  is  the  value  of  the  state 
which  maximizes  the  probability  P(r).  Since  the  error  statistics  are  assumed  to  be  Gaussian 

P(r)  =  (271)-'"'^  |Q|"'^  exp|^-|r"Q-'r  j  (26) 

where  the  residual  vector  r  is  the  residual  vector  associated  with  the  estimate  x .  Substituting  Eq  (21) 
into  Eq  (26)  leads  to  the  minimization  of 
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0r-'ISi)'Q-'(r-T5i)]  =  O 


where  r  is  the  total  residual  vector  associated  with  the  reference  state  From  here  it  is  pos^e  to 
solve  for  the  correction,  new  estimate  and  associated  covahaiKX  matrix  as 

5x(/o)  =  (T^Q'T)'‘T^Q-‘r  (28) 

i(/o)  =  x„,(/o)  +  5x(/o)  (29) 

^81  “  ^(*8i*8i  ) 

/  x-i 

p,  =(rQ-'T) 

where  T,  Q,  and  r  are  the  total  matrices  and  vector  defined  earlier  (4:60-63,70). 

It  is  important  to  note  that  the  estimate  found  in  Eq  (29)  is  the  correct  answer  only  if  the  reference 
has  converged  to  a  solution.  In  general,  the  estimate  becomes  a  new  reference  state  and  the  whole  process 
is  repeated  until  convergence  is  attained.  One  practical  convergence  standard  is  for  successive  values  of 
5x  to  be  within  the  one  o  values  provided  by  P5x  and  for  the  residuals  to  also  agree  with  their  expected  a 

values  (4:  77). 

2.6  Implementing  Least  Squares  on  a  Computer 

When  working  with  large  amounts  of  data,  the  matrix  and  vector  equations  introduced  above  can 
become  large  and  unwieldy.  For  this  reason,  it  is  usual  to  process  the  data  one  observation  at  a  time  using 


the  following  miming  sums: 
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After  all  of  the  data  is  processed  the  correction  to  the  state  and  the  covariance  matrix  can  be  fiHmd  from 


Pfc  =  A-' 

5*(fo)  =  Pa.B 


(33) 


where  Sx  is  the  correction  that  turns  the  referance  state  into  either  the  estimate  or  a  better  reference  state 
and  Psx  is  the  covariance  matrix  associated  with  the  corrected  state  (4:64,72). 


27  Summary 

The  method  of  least  squares  is  a  proven  approach  to  estimating  the  state  of  a  satellite  from  a  large 
mimber  of  observations.  When  used  in  a  differential  corrector  it  allows  for  the  iterative  correction  of  an 
initial  reference  state  to  the  correct  estimate  of  the  true  state.  It  does  this  by  adjusting  the  reference  state 
to  produce  a  tnyectory  that  agrees  best  (i.e.  minimizes  the  squared  error)  with  all  of  the  data  instead  of 
agreeing  perfectly  with  only  two  or  three  observations.  Least  squares  also  has  the  property  that  the 
accuracy  of  the  estinutte  increases  as  the  number  of  data  points  increases  although  the  relationship  is  not 
linear.  For  example,  if  only  one  variable  is  involved  it  can  be  shown  that; 


where  the  standard  deviation  of  the  estimate,  N  is  the  number  of  data  points,  and  oj  is  the  standard 
deviation  of  each  measurement  (a  constant  in  this  case). 
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III.  Methodology 


3.1  Introduction 

The  data  used  in  this  research  was  produced  using  a  two-step  procedure  in  which  simulated 
observation  data  was  converted  into  a  state  vector  and  associated  covariance  matrix.  The  first  step 
involved  the  use  of  a  truth  model  to  generate  simulated  range,  azimuth,  and  elevation  observatkHis  for  a 
single  track  of  a  satellite  as  seen  by  a  space  surveillance  radar  located  on  the  surfsce  of  the  earth.  In  the 
second  step,  this  track  of  data  was  reduced  to  a  state  vector  and  covariance  matrix  using  least  squares 
differential  correction.  Error  statistics  were  also  generated  and  the  accuracy  of  the  estinuited  state  was 
calculated  by  comparing  it  to  the  true  state. 

Analysis  of  the  data  concentrated  on  evaluating  how  good  the  estimator  was  in  concentrating  all  of 
the  data  available  in  the  track  into  the  estimated  state  vector.  The  effects  of  different  dynamics  models, 
data  acquisition  rate,  and  maximum  elevation  of  the  pass  were  investigated  using  the  error  statistics  and 
accuracy  of  each  estimate. 

3.2  The  Truth  Model 

The  truth  model  simulates  the  orbital  motion  of  a  satellite  by  using  two-body  dynamics  pixis  J2  2onal 
harmonic  and  atmo^heric  drag  perturbations.  Because  this  study  concentrates  on  relatively  short  arcs  (22 
minutes  or  less)  of  near  earth  orbits,  it  was  felt  that  other  perturbations  such  as  the  higher  order 
geopotential  terms  and  n-body  effects  were  insignificant  and  could  be  ignored. 

3. 2. 1  Program  Execution.  Input  to  the  truth  model  consists  of  an  initial  state  containing  position  and 
velocity  (expressed  in  ECI  coordinates),  an  epoch  (or  start)  time,  the  end  time  for  the  arc  in  question,  and 
the  time  step  between  each  data  point.  The  program  then  numerically  integrates  the  equations  of  motion 
using  a  fourth  order  Hamming  predictor-corrector  routine.  The  predictor  part  is  based  on  a  fourth-order 
Milne  method  predictor  and  the  corrector  method  is  derived  from  fourth-order  Milne  and  Adams-Moulton 
predictor-corrector  methods  (5:391).  Since  this  method  requires  four  previous  values  of  the  state  in  order 
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to  predict  the  next  value,  a  Picard  iteration  is  used  to  generate  three  future  state  values  based  (m  the  initial 
state  (6:108-109). 

OuQNit  generated  by  the  program  consists  of  qrhemerides  and  observations  which  are  written  to  two 
separate  files  for  later  use  by  the  estimator  programs.  After  each  ephemeris  poim  is  generated  a 
subroutine  is  called  to  convert  the  ECl  position  of  the  satellite  into  topocentric  observations  relative  to  an 
imaginary  sensor  site  located  in  the  Indian  Ocean  (see  A|^ndix  A).  If  the  elevation  of  this  "observation" 
is  greater  than  zero  then  it  is  assumed  that  the  sensor  can  see  the  satellite  and  so  the  ephemeride  and  its 
associated  observation  are  written  to  their  reflective  data  files.  If  the  elevation  is  less  than  zero,  then  the 
ouqjut  data  files  are  not  modified  and  the  next  position  calculated.  This  basic  cycle  is  repeated  until  the 
specified  end  time  has  been  reached. 

In  order  to  better  model  the  real  world,  the  observations  are  tiKxlified  by  the  addition  of  simulated 
errors  or  noise.  This  is  done  using  a  Gaussian  pseudo-random  number  generator  and  assumed  values  for 
the  one-sigma  accuracy  of  the  radar  in  range,  azimuth,  and  elevation.  For  purposes  of  this  stiufy  the 
assumed  values  were;  1(X)  meters  in  range,  0.02S  degrees  in  azimuth,  and  0.025  degrees  in  elevation. 

Since  a  new  seed  value  is  used  every  time  the  number  generator  is  called,  each  value  output  will  be 
different  and  the  overall  distribution  will  be  very  close  to  the  normal  distribution  shown  in  Figure  1. 
Multiplying  the  one-sigma  values  by  the  random  numbers  results  in  the  simulated  errors  also  following  a 
Gaussian  distribution. 

Several  special  options  were  built  into  the  truth  model  program  to  facilitate  the  testing  and  research 
objectives.  The  program  allows  for  the  time  step  to  be  changed  from  within  the  program.  This  meant  that 
the  input  file  did  not  have  to  be  edited  for  the  single  imrpose  of  changing  the  time  stq).  Also,  the  addition 
of  noise  or  drag  perturbations  ctHild  be  turned  off  prior  to  starting  the  numerical  integration.  Turning  off 
the  noise  permits  the  generation  of  "perfect"  data  for  testing  purposes,  while  turning  off  the  drag  allows 
for  investigation  of  the  significance  of  drag  in  affecting  the  estimate  and  covariance  matrix. 
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3.2.2  The  Equations  of  Motion.  The  baseline  for  the  equations  of  motion  is  the  first  order  difTerential 
equations  from  the  solution  of  the  two-bo(fy  problem.  In  vector  notation,  the  two-body  equatitm  of  motion 

is: 

r  +  -^r  =  0  (34) 

r 

where  the  vector  r  extends  from  the  earth':,  'xnter  of  mass  to  the  satellite's  cemer  of  mass,  r  is  the 
magnitude  of  this  vector  and  p.  is  the  gravitational  parameter  of  the  earth  and  is  defined  as  psGA/  , 
where  G  is  the  universal  gravitational  constant  and  A/ is  the  mass  of  the  earth.  In  a  rectangular  coordinate 
system  (such  as  ECI),  Eq  (34)  leads  to  six  first  ortter  differential  equations; 

r  =  (35) 

x  =  u 

y=^v  (36) 

z  =  w 

u  =  -\ixfr^ 

v  =  -\iylr^  (37) 

w  =  -\u/r^ 

where  x.y.zsae  the  components  of  the  position  vector,  and  ».  v,  w  are  the  components  of  the  velocity 
vector  (reference  Figure  2). 


Figure  2.  Position  and  velocity  components  in  ECI  coordinates. 
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The  most  significant  of  the  various  geopotential  efiSxts  is  found  in  the  second  zonal  hannonic,  usually 
referred  to  as  the  J2  getqwtential  term  or  simply  the  J2  term,  which  is  caused  the  earth's  oUateness 
(6:59,80).  This  term  is  normally  written  in  qjherical  coordinates  as: 

r„=i!^(3cos=e-i)  (38) 


where  is  the  earth's  equatorial  radius,  ^2  is  coefficient  for  the  second  zonal  harmonic,  and  0  is  the 
angle  between  the  Z  coordinate  axis  and  the  satellites  position  vector  (reference  Figure  3)  (6:59).  In  order 
to  convert  this  geopotential  term  to  ECI  coordinates  it  is  necessary  to  rewrite  the  angle  as  cos6  =  z/r.  If 
the  constants  on  the  right  side  of  Eq  (38)  are  collected  as: 

^  =  (39) 

then  Eq  (38)  can  be  rewritten  in  cartesian  coordinates  as: 


20 


(40) 


z 


Figure  3.  Coordinates  for  Geopotential 


The  effect  of  the  geopotential  manifests  itself  as  an  additional  acceleration  of  the  satellite  which  is 
related  to  the  geopotential  through  the  gradient  operator: 

*20  (41) 
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where  *20  ^  ^  additional  acceleration  due  to  the  J2  geopotential  term.  Calculating  the  gradient  results 
in  an  additional  term  for  each  the  velocity>dot  equations; 


*  ~  ■*‘*20 


'3x/r*-15xr*  /r"' 

V 

y 

-A 

-\Syz^  tr'' 

r 

^  9zlr^ -\Sz^  Ir''  ^ 

The  basic  equation  for  modeling  the  effect  of  atmoq^ric  drag  is; 


2  m 

|=|{i+e>'.>'-ex,i}| 


(43) 


where  a^  is  the  effect  of  air  drag,  Q)  is  the  satellite's  coefficient  of  drag.  A/m  is  the  satellite's  ratio  of 
cross  sectional  area  to  mass,  p  is  the  atmo^heric  density,  is  the  satellite's  velocity  relative  to  the 

atmosphere,  and  6  is  the  rotational  rate  of  the  earth  (atmo^here  is  assumed  to  rotate  with  the  earth) 
(2:423-424).  The  quantity  C^A/m  is  often  lumped  together  and  called  5*  .  For  the  purposes  td’ this 
shuty,  the  satellite  was  assumed  to  have  a  C£)  =  2.Q,A  =  7.5  m^,  and  m  =  1000  kg. 

The  atmospheric  density  was  modeled  using  the  1976  Standard  Atmoq>here  in  which  the  density,  as  a 
ftinction  of  height,  is  given  by: 


p(/i)  =  exp 


J flV  +( flV  +( — T 

v/tj  1^9/ 


where  p(h)  is  in  kg-m’^,  A  is  the  altitude  in  kilometers  and 

zi  =  246482.873  22  =  5377132.030  23  =  1536228.60  24  =  19174788.0  25  =  105.381650 

26  =  4%8.56120  27  =  56060.0780  2g  =  160314.760  29  =  344944.780 

This  nmdel  is  valid  for  altitudes  between  100-1000  km  (8). 
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3.2.3  Program  Validation.  Proper  operation  of  tlK  truth  OKxkl  was  validated  by  inciemratally 
testing  its  accuracy  using  three  models  of  the  dynamics;  (1)  two-body  only,  (2)  two-body  plus  J2,  and  (3) 
two-body  plus  J2  plus  drag.  The  two-body  test  was  done  using  a  circular  reference  orbit  whidi  allowed  for 
easy  calculation  of  the  satellites  position  as  a  function  of  time.  To  test  the  J2  part,  two  different  reference 
orbits  were  used  (see  Table  1). 


Table  1 .  Reference  Orbits  for  Truth  Model  Validation 


Inclination 

(deg) 

Eccentricity 

Period 

(min) 

Apogee  Ht 
(km) 

Perigee  Ht 
(km) 

Nodal 

Regression 

Apsidal 

Rotation 

7.33 

0.0166 

91.19 

444.4 

221.6 

8.3  deg/day 

63.4 

0.05 

103.43 

1285.5 

555.6 

Each  orbit  was  propagated  for  four  days  in  order  to  observe  the  changes  in  the  right  ascension  of  the 
ascending  node  and  the  argument  of  perigee.  Since  the  rates  shown  in  Table  1  agreed  very  closely  with 
the  predicted  values  it  was  felt  that  the  J2  equations  were  validated  (2;  1S7-1S8). 

3.3  The  Estimator 

The  estimator  is  a  least  squares  differential  corrector  using  a  dynamics  model  based  on  two-body  and 
J2  effects.  By  using  a  Taylor  Series  afvroximation,  the  equations  of  state  and  the  state  transition  nurtrix 
are  expressed  as  explicit  functions  of  time.  Four  different  dynamics  models  are  investigated  by  varying 
the  order  of  the  approximations  used  for  the  two-bo  ly  a^d  the  J2  parts  of  the  equations.  These  variations 
are  detailed  in  Table  2. 


Table  2.  Highest  Order  of  Approximation  in  Equations  of  State 


Estimator 

Two-body 

Position  /  Velocity 

h 

Position  /  Velocity 

ets4 

4th  /  3rd 

None 

ets4a 

4th /3rd 

2nd /1st 

ets4b 

4th /3rd 

2nd  +  3rd  /  1st  +  2nd 

ets5 

5th /4th 

2nd  +  3rd  /  1st  +  2nd 

A  different  state  transition  matrix  is  derived  from  each  set  of  state  equations,  except  for  'etsS'  which 
uses  the  same  state  transition  matrix  as  'ets4b'. 
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3. 3. 1  Program  Execution.  The  input  to  the  estimator  program  consists  of  the  observation  and 
qthemeris  files  generated  by  the  truth  program.  Before  the  estimation  process  actually  begins,  all  the 
observations  are  read  into  an  array  for  later  use.  During  this  process  the  program  determines  how  many 
observations  are  in  the  file.  By  storing  all  of  the  observations  in  an  array,  the  program  does  not  need  to 
read  the  observation  file  more  than  once  thus  reducing  program  run  time.  This  approach  also  allows  fm 
the  estimate  to  be  epoched  to  the  time  of  the  middle  observation  which  effectively  doubles  the  time  span 
over  which  the  Taylor  Series  approximations  to  the  equations  of  state  are  effective.  After  determining 
which  observation  is  the  epoch  observation,  the  corresponding  ephemeride  can  be  lead  thus  providing  the 
true  state  at  epoch. 

A  very  simple  technique  is  used  to  determine  the  initial  reference  state  at  the  beginning  of  the 
estimation  process.  The  position  part  of  the  reference  state  is  immediately  available  from  the  qxxh 
observation  by  converting  the  topocentric  values  to  ECI  coordinates  (see  Aiq)endix  A).  The  velocity  part 
of  the  reference  state  is  then  found  using  two  consecutive  observations  (epoch  and  epoch+l)  by  dividing 
the  change  in  each  position  component  the  time  between  observations; 


v  = 


X  _ X 

tpoch 

^»fodn\  ^tpoeh 


^*pocln-\  ^tpoch 
^»poch+l  ~  ^tpoch 
^*poch*\  ^tpoch 


(44) 


The  estimator  uses  an  iterative  aiqiroach  to  correct  the  initial  reference  state  into  an  estimate  of  the 
satellite's  state  at  the  epoch  time.  At  the  beginning  of  each  iteration,  the  program  calculates  several 
"constants”  based  on  the  reference  state  for  use  in  the  equations  of  state  and  state  transition  matrix 
subroutines  and  initializes  the  running  sum  matrix  and  vector  defined  in  Eq  (3 1)  and  Eq  (32). 
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As  each  observation  is  processed,  a  predicted  observation  is  obtained  by  evaluating  the  equations  of 
state  using  the  difference  between  the  observation  time  and  the  epoch  time  ( dt  =  •  tg^^ ) 

and  then  converting  from  ECI  to  topocentric  coordinates.  The  observation  is  then  compared  with  the 
predicted  observation  and  the  residual  calculated.  The  residuals  for  range,  azimuth,  and  elevation  are  also 
used  to  tqxiate  ruiming  sums  for  the  mean  and  root  mean  squared  (RMS)  residuals; 

N 

sum  =  '^residua^ 
sum2  =  ^{residua^ ) 

i=l 

The  predicted  state  is  next  used  to  calculate  the  data  linearization  matrix  H  (see  ^)pendix  B)  which 
is  then  combined  with  the  state  transition  matrix  to  produce  the  observation  linearization  matrix  T  =  H<I>. 
After  multiplying  by  the  instrumental  covariance  matrix  Q  and  the  residual  vector  r,  the  results  ate  added 
to  the  running  matrix  and  vector  sums  found  in  Eqs  (3 1)  and  (32). 

When  all  of  the  observations  have  been  processed,  a  new  estimate  of  the  state  and  its  covariance 
matrix  P  is  calculated  using  Eqs  (29)  and  (33).  Also  the  mean  and  RMS  residuals  ate  calculated  for  the 
previous  estimate: 


mean  residual 
RMS  residual 


sum 
obcnt 
sum! 
V  ohcnt 


(46) 


where  obcnt  is  the  total  number  of  observations,  sum  is  the  cumulative  residuals  for  all  observations,  and 
sum2  is  the  cumulative  residuals  squared  for  all  observations. 

If  the  change  in  each  element  of  the  state  is  smaller  than  one-^Krcent  of  its  corresponding  diagonal 
entry  in  the  covariance  matrix,  then  the  estimator  is  considered  to  have  converged.  If  not,  then  the 
estimate  becomes  a  new  reference  state  and  the  program  loops  back  to  start  another  iteration. 
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After  convergence,  the  program  calculates  a  new  set  of  residual  statistics  by  compariiig  each 
observation  to  predictions  based  on  the  new  estimated  state.  The  accuracy  of  the  estimate  is  also 
calculated  using  the  final  estimated  state  and  the  true  state  to  find  the  magnitude  oi  the  position  and 
velocity  errors  at  qx>ch: 

)'  +  y 

\ -  (47) 

3.3.2  Development  of  the  Taylor  Series  Approximations.  As  mentioned  earlier,  both  the  equations  of 
state  and  the  state  transition  matrix  are  approximated  using  a  Taylor  Series  expansion  in  time.  The 
primary  advantage  of  taking  this  a|q>roach  rather  than  numerically  integrating  the  equations  (tf  motion 
and  equations  of  variation  is  that  the  state  and  state  transition  matrix  can  easily  be  calculated  at  any  given 
time.  It  was  also  felt  that  working  with  explicit  fimctions  of  time  would  reduce  the  computational  burden 
placed  on  the  computer  -  which  is  one  of  the  caveats  mentioned  in  Chapter  1 . 

To  start  with,  the  Taylor  Series  rqrresentation  of  the  position  and  velocity  of  a  satellite  can  be  written: 


W  0  0  2  0  6  °  24  »  120  » 

v(/)  =  Vq  ■^7*o<**  +0(t**) 

2  O  24 

where  r(t)  is  the  satellite's  position  at  some  time  t,  v(t)  is  the  satellite's  velocity  at  some  time  t,  tq  is  the 
satellite's  position  at  t  =  0 ,  vq  is  the  satellite's  velocity  at  t  -  0 ,  and  sq  is  the  satellite's  acceleration  at  t  = 

0 .  The  notation  of  using  a  subscript  zero  to  indicate  a  condition  at  t  =  0  will  trow  be  dropped  and  it  will 
be  understood  that  all  quantities  are  evaluated  at  t  =  0 .  As  was  done  in  section  3.2.2,  the  acceleration  of 
the  satellite  can  be  brcricen  into  two  parts,  the  two-body  term  and  the  J2  term: 

*  ~  ^2body  ■*■•20  (49) 


Rewriting  Eq  (34)  as 


24 


(50) 


produces  the  second  order  term  for  the  two-body  dynamics.  Differentiating  with  respect  to  time  gives 


a 


ujr^  -3x(«3c  +  ty +  H'z)/r’  '' 
vlr^  - 3y(ux  +  yy  +  wz)/r* 
w/r^  -3z(ux  +  yy +  wz)/r^  ^ 


(51) 


which  is  the  third  order  two-body  term. 


a 


lr^-3x(y^  +rda)lr^ -6u{rdv) I -^-XSxirdvf 
/  -  3y(V^  +  rda)  / r*  -  6v{rdv)  /  r*  + 1  Sy{rdv)^  I  r'' 

/  -  3z(V^  +  rda)  / r’  -  6w{rdv)  / r*  + 1  Szirdvf  lr\ 


(52) 


which  is  the  fourth  order  two^xxly  term  where 


rdv  =  ux-\rvy  ■¥y^z 
rda  =  xa^  +)Kiy  +za^ 

V  s 


and  finally 


a 


( 


V 


-Uu/r‘  +45[M(rdv)’  +*(rrfv)(F’  -4/r)]/r’  -9«(F’  -n/r)/r’  +  15|ix(n/v')/r'  -lOSxCrrfv)’  /r*' 
-Uv/r*  +45[v(n#v)’  +>'(«/v)(F"  -n/r)]/r’  -9v(F’  -p/r)/r’  +  15p>'(rrfv)/r*  -105y(«/v)’  /r 
-y^lr*  ^ 45{ w(rrfi;)’  +  r(r<#w)(F’  -  m  / r)] / r’  -  9w(F'  -  |i  / r)  / r’  + 1 5ur(nfr) / r*  -  1052(«fr)’  / r 


(53) 


which  is  the  fifth  order  two-body  term  where 
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(54) 


a,  =  -|ix/r^ 
a^=-^ylr^ 
a,  =  -)u/r* 

Refisrring  back  to  Eqs  (39-42),  the  perturbative  acceleration  due  to  J2  effects  can  be  written  as; 


a 
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=  -y4 


- \Sxz^ 

“iy/r^  -15^^  /r’ 

^  9r/r*-15z'/r’  ^ 


udiich  is  the  first  order  J2  term  for  the  Taylor  Series. 


(55) 


Differentiating  with  req)ect  to  time  gives; 


'3tt/r*-15 

2hocz  +  +  x(rdv) 

/r’+105xz 

1 

II 

8 

3v/r’  -15 

iMyz  +vz^  +y(rdv) 

/r’+105>«*(r</i')/r’ 

^  9w/r* 

-  45[wz *  - zirdv)] Ir''  105r' (rdv)  / r’  J 

which  is  the  second  order  Taylor  Series  term  for  the  J2  effect. 

Develq>ment  of  the  state  transition  matrix  equations  follows  directly  fiom  the  equations  of  state  as; 

=  (14) 

where  the  fimction  h  is  given  by  Eq  (48)  as  an  explicit  fiinction  of  time.  A  complete  derivatitm  of  the 
state  transition  matrix  is  given  in  Appendix  C.  It  should  be  noted  that  the  estimate  is  much  more  sensitive 
to  the  order  of  the  equations  of  state  than  it  is  to  the  order  of  the  state  transition  matrix  (9;  161). 

3.3.3  Program  Validation.  The  estimator  program  was  validated  by  using  it  to  estimate  the  state  for 
orbital  arcs  taken  fiom  several  test  orbits.  For  these  tests,  the  addition  of  noise  and  drag  perturbations  was 
"turned  off"  so  that  proper  operation  of  the  program  could  be  verified.  The  orbital  parameters  for  the  six 
test  orbits  used  are  given  in  Table  3. 
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Table  3.  Test  Ort>its  fw  Estimator  Validation 


Orbit  ID 

Eccentricity 

A 

35.26 

0.3772 

412.4 

18895 

5053 

B 

4.33 

0.0010 

89.45 

254 

241 

C 

93.06 

0.0010 

89.45 

254 

241 

D 

35.26 

0.0819 

98.06 

1243 

89 

E 

0.00 

0.0000 

88.93 

222 

222 

F 

7.33 

0.0166 

91.19 

444 

222 

The  major  program  nnodules  that  make  up  the  program  were  also  individually  tested.  This  included  the 
equations  of  state,  the  state  transition  matrix,  the  time  conversion  routines,  the  data  linearization  matrix 
H,  and  the  routines  to  convert  from  topocentric  to  ECI  and  back. 

The  equations  of  state  were  tested  by  computing  the  state  of  a  satellite  at  three  second  intervals  using 
the  Taylor  Series  equations  and  then  comparing  with  the  numerically  integrated  state  using  the  equations 
of  motion  developed  for  the  truth  model.  The  state  transition  matrix  was  also  tested  comparison  to  the 
numerically  integrated  matrix  but  only  at  two  times  (dt  -  3  sec  and  dt  =  S7  sec).  The  matrix  was  first 
tested  based  only  on  two-body  dynamics  and  once  that  was  validated  the  Jj  terms  were  added.  In  order  to 
test  the  J2  terms,  the  value  of  the  J2  coefficient  was  temporarily  changed  to  be  the  same  as  the  earth's 
gravitational  parameter  so  that  the  effects  could  more  easily  be  seen.  In  order  to  get  good  agreement 
between  the  approximate  and  "exact"  matrix  it  was  necessary  to  use  a  third  order  expression  for  the  J2 
perturbation.  Of  course,  this  degree  of  accuracy  was  not  really  needed  since  the  true  value  of  the  J2 
coefficient  is  much  smaller  than  the  value  for  |x. 

3.3.4  Determining  the  Order  of  the  State  Equations.  Initially  the  two-body  dynamics  were  modeled 
using  a  third  order  and  a  fourth  order  set  of  equations  to  see  if  there  was  much  difference  between  them. 
The  tests  shown  in  Table  4  demonstrated  that  the  fourth  order  equation  was  significantly  better  and  so  it 
was  used  as  the  basis  for  the  equations  of  state  and  state  transition  matrix.  For  this  test,  the  Taylor  Series 
equations  included  only  two-body  effects  while  the  truth  data  was  generated  using  a  numerical  integrator 
with  two-body,  J2  and  drag  effects.  Errors  were  calculated  by  computing  the  vector  magnitude  of  the 
difference  between  the  two  states. 
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Table  4.  Position  and  Velocity  Errors  from  3rd  and  4th  Older  Taylor  Series 


Test 

3rd  Older  T.S. 

WKSSS^BM 

# 

Velocity  Error  (cm/sec) 

1 

0.092 

0.381 

0.007 

0.163 

2 

0.101 

0.431 

0.004 

0.130 

3 

0.143 

0.947 

0.004 

4 

0.107 

0.968 

0.002 

0.100 

5 

0.025 

0.589 

0.0001 

0.017 

a.  Tests  1  and  2  were  run  using  test  orbit  E  from  Table  3. 

b.  Tests  3-S  were  run  using  test  orbit  F  from  Table  3. 


After  developing  the  first  and  second  order  equations  for  the  J2  perturbation,  a  different  test  was  run 
to  verify  the  equations  and  to  provide  some  insight  into  their  relative  performance  differences.  For  this 
test,  a  reference  orbit  was  propagated  using  the  Hamming  numerical  integrator  and  a  time  stq)  of  three 
seconds.  This  "true"  value  of  the  state  was  then  compared  to  the  predictions  of  three  different  equations  of 
state  and  the  results  output  as  the  magnitudes  of  the  position  and  velocity  errors  (Rmag  and  Vmag).  The 
three  different  equations  of  state  were;  (1)  The  fourth  order  two-body  equations  only  ('eom4.for').  (2)  The 
fourth  order  two-body  equations  plus  the  first  order  J2  equations  ('eom4a.for*).  (3)  The  fourth  order  two- 
body  equations  plus  the  first  and  second  order  J2  equations  ('eom4b.for').  When  the  truth  model  included 
only  the  two-body  dynamics  'eom4'  performed  the  best.  However,  with  the  addition  of  J2  and  drag  effects, 
it  was  quite  obvious  that  the  S2  equations  would  be  needed  in  the  estimator.  Figure  4  shows  the  relative 
performance  of  the  three  equations  of  state  as  a  plot  of  their  position  error  (Rmag)  versus  time  (in 
timesteps). 
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Figure  4.  Position  Error  vs.  Timestq)  for  eom4,  eoin4a,  and  eoin4b. 

The  data  for  Figure  4  was  generated  using  test  orbit  A  from  Table  3.  As  can  be  seen,  the  plots  show 
that  a  significant  performance  improvement  occurs  when  J2  is  taken  into  account.  There  is  very  little 
difference  between  'eom4a'  and  ‘eom4b'  but  'eom4b'  is  slightly  better.  This  is  expected  since  'eotn4b'  uses 
a  higher  order  approximation  for  J2  effects. 

However,  for  very  low  orbits,  such  as  test  orbits  B  through  F  in  Tsmle  3,  'eom4a'  actually  performs 
better  than  'eom4b'.  As  can  be  seen  in  Figure  5  (generated  from  orbit  C  in  Table  3),  'eom4b'  is  initially 
more  accurate  than  'eom4a',  but  after  about  100  seconds  'eom4a'  actually  has  a  smaller  error.  This  sort  of 
behaviour  is  unexpected  since  it  is  normally  expected  that  a  higher  order  of  approximatirm  would  give 
more  accurate  results. 
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Figure  3.  Position  Error  vs.  Timestep  for  eom4a  and  eoin4b. 

Using  a  figure  of  merit  that  the  position  error  should  be  within  one  kilometer,  the  duration  during 
which  each  version  of  the  state  equation  should  be  usable  is  given  in  Table  S. 


Table  5.  Effective  Time  Duration  for  State  Equations 


Orbit 

eom4 

eom4a 

eom4b 

E 

5:33 

6:12 

6:00 

A 

9:18 

9:21 

B 

5:48 

6:36 

6:24  1 

As  can  be  seen,  for  very  low  orbits,  'eom4a'  is  slightly  better  than  ’eom4b'  and  both  of  them  are 
significantly  better  than  'eom4'.  However,  at  higher  altitudes  the  difference  between  the  three  versions  is 
very  slight  and  the  an>roximate  equations  of  state  are  "good"  for  50-70%  longer. 
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3.4  Orbital  Parameters. 


All  of  the  orbits  mentioned  so  far  were  used  only  for  develc^ment  and  testing  of  the  programs  used  in 
this  study.  Analysis  of  the  estimator  performance  was  actually  carried  out  using  the  oibits  detailed  in 
Tables  6  and  7. 


Table  6.  Orbital  Elements  -  Satellite  20003 


1  ECl  State  Vector 

Classical  Elements  I 

Epoch 

Variable 

Period 

95 

xpos 

4%1.174 

km 

Eccentricity 

0.0001 

ypos 

-4210.369 

km 

inclination 

28.5 

deg 

zpos 

-2286.044 

km 

Rt.  Ascension 

0 

deg 

X  vel 

5.280 

km/sec 

56 

deg 

y  vel 

km/sec 

True  Anomaly 

260 

deg 

z  vel 

2.610 

km/sec 

Perigee  Ht 

517.9 

km 

519.3 

km 

Table  7.  Orbital  Elements  •  Satellite  20004 


1  ECI  State  Vector 

Classical  Elements  I 

Epoch 

Variable 

Period 

no 

min 

xpos 

1968.306 

km 

Eccentricity 

0.0001 

ypos 

-6455.632 

km 

Inclination 

28.5 

deg 

zpos 

-3505.122 

km 

Rt.  Ascension 

0 

deg 

X  vel 

6.993 

km/sec 

25 

deg 

y  vel 

1.647 

km/sec 

True  Anomaly 

260 

deg 

zvel 

0.894 

km/sec 

1225.9 

km 

1227.4 

km 

The  very  low  eccentricity  is  typical  for  low  earth  satellites  that  have  been  in  orbit  for  any  significant 
time  and  also  eliminates  the  influence  of  apogee  and  perigee  on  the  length  of  each  pass.  The  inclination 
was  chosen  because  it  is  typical  of  U.S.  shuttle  orbits  and  because  it  was  high  enough  to  ensure  a  ninety 
degree  pass  over  the  sensor  site  used  in  the  study.  The  argument  of  perigee  was  chosen  to  put  the  perigee 
in  the  northern  hemisphere  which  is  where  the  sensor  site  is  located.  By  setting  the  true  anomaly  to  260° 
the  satellite  is  positioned  for  an  ascending  pass  and  will  begin  its  trajectory  outside  the  sensor's  coverage. 
This  reduces  the  length  of  the  orbit  which  must  be  integrated  before  the  satellite  enters  sensor  coverage. 
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Satellite  20003  was  used  to  evaluate  the  estimator  performance  on  very  low  ME  satellites.  Track 
duration  varied  from  about  6.5  minutes  at  very  low  elevation  to  12.4  minutes  at  very  high  elevatkm. 
Satellite  20004  was  chosen  to  evaluate  estimator  performance  on  relatively  high  NE  satellites.  Track 
duration  varied  from  10.2  minutes  to  2 1 .5  minutes.  The  inclination  and  eccentricity  were  kqit  constant  so 
that  the  effects  of  altitude  could  be  more  clearly  noted. 

3.5  Parameters  for  Data  Generation. 

Three  major  parameters  were  used  in  structuring  the  data  generation;  (1)  the  duration  of  the  pass,  (2) 
the  data  rate  in  observations  per  second,  and  (3)  the  dynamics  used  in  the  estimator.  The  maximum 
elevation  of  a  track  is  directly  related  to  how  long  the  satellite  is  visible  to  the  site  and  therefore  the 
amount  of  data  available  for  the  estimator.  This  relationship  is  not  linear  but  follows  the  pattern  shown  in 
Figured. 


For  a  given  sensor  site,  the  exact  relationship  between  the  length  of  pass  and  maximum  elevation  is 
dependent  on  the  satellite's  orbit  (shape  and  size)  but  will  still  have  the  basic  shape  exhibited  in  Figure  6. 
In  order  to  characterize  the  relationship  between  estimator  performance  and  pass  length,  data  was 
collected  using  the  schedule  shown  in  Table  8.  As  shown  in  Figure  6,  the  duration  of  a  pass  changes  more 
quickly  at  low  elevation  than  at  high  elevation  and  after  about  30  degrees  there  is  relatively  very  little 
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change  at  all.  Because  of  this,  more  passes  were  scheduled  at  the  lower  elevations  than  at  the  higher 
elevatitms. 


Table  8.  Schedule  of  Passes  for  Data  Collection 


Orbit 

Para 

m 

H 

B 

m 

D 

E 

H 

G 

H 

1 

J 

20003 

Time 

6.60 

7.56 

8.37 

9.65 

10.41 

11.36 

12.08 

12.28 

12.45 

n/a 

M 

El 

3.45 

4.84 

6.35 

EC3 

13.16 

20.47 

35.38 

48.03 

89.44 

n/a 

20004 

Time 

10.2 

12.85 

15.05 

16.35 

17.50 

18.90 

19.95 

20.95 

21.30 

21.50 

n 

El 

3.60 

6.28 

9.77 

12.64 

16.03 

22.36 

KiT!!! 

47.24 

65.22 

87.72 

a.  Time  refers  to  the  duration  of  the  scheduled  pass  in  minutes. 

b.  El  refers  to  the  max  elevation  of  the  scheduled  pass  in  degrees. 


The  rate  at  which  a  space  surveillance  radar  generates  observations  is  influenced  by  several  fectors: 
the  type  of  radar,  the  radar's  operating  frequency,  and  the  range  to  the  satellite.  This  makes  it  difflcult  to 
specify  a  single  data  rate  as  being  representative  of  all  the  radars  in  the  SSN.  However,  it  can  be  assumed 
that  a  radar  could  easily  be  capable  of  generating  at  least  two  observations  per  second  and  perhaps  as 
many  as  100  observations  per  second  (7.156).  It  should  be  noted  that  the  maximum  data  rate  that  sensors 
are  expected  to  transmit  to  the  SSC  is  one  observation  every  six  seconds  (10).  Based  on  this  information, 
data  for  this  study  was  generated  at  two  rates  -  one  observation  every  six  seconds  and  two  observations 
per  second. 

3.6  Experimental  Procedure. 

The  first  step  was  to  determine  the  proper  epoch  from  which  to  start  the  satellite's  trajectoiy  in  order 
to  obtain  the  desired  maximum  elevation  during  each  pass.  This  was  done  using  a  program  written 
expressly  for  this  purpose.  The  output  from  this  program  gives;  (1)  an  epoch  for  the  state  vector  in  Table 
5,  (2)  the  associated  rise,  culmination,  and  set  times  and  elevation  angles,  and  (3)  the  satellite's  state  and 
epoch  one  minute  before  entering  sensor  coverage. 

Next,  the  resulting  state  vector  and  epoch  was  input  into  the  truth  model.  By  using  the  set  time  plus  a 
safety  margin  the  end  time  for  the  truth  model  propagation  could  be  efficiently  specified.  The  time  step 
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was  q)ecified  either  in  the  iiqHit  file  or  during  program  execution.  The  ou^t  from  the  truth  model  was 
then  iqxit  into  the  four  different  estimators  for  processing. 

The  data  generated  by  the  estimators  was  then  plotted  to  allow  analysis  and  comparison  of  each 
estimator's  performance. 

3.7  Equipment. 

The  programs  used  in  this  study  were  written  in  Microsoft  FORTRAN  vS.l  and  run  on  a  486DX  PC 
running  at  33  MHz.  Some  supplementary  work  and  generation  of  graphs  was  done  using  Mathcad  v3. 1. 
Mathematica  was  used  extensively  in  deriving  the  Taylor  Series  equations  for  the  equations  of  state  and 
the  state  transition  matrix. 

3.8  Summary. 

The  effectiveness  of  using  a  least  squares  batch  estimator  to  estimate  the  state  of  a  satellite  is 
investigated  using  simulated  observations  from  a  single  radar  tracks.  The  observations  are  generated 
using  a  truth  model  that  uses  a  fourth  order  numerical  integrator  to  prqragate  the  satellite's  orbit  for  a 
small  portion  of  the  orbit.  The  estimator  models  the  orbital  dynamics  using  only  two4xxly  and  J2  effects 
using  equations  of  state  and  state  transition  matrix  that  are  expressed  in  explicit  form  using  a  Taylor 
Series  rqrproximation. 
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IV.  Results  and  Analysis 


4.1  Introduction 

Four  different  versions  (ff  the  estimator  were  compared  using  data  from  two  different  orbits,  two 
different  time  steps,  and  two  different  noise  settings  (on  and  off).  The  results  were  examined  to  determine 
whether  (1)  the  estimators  were  effectively  conqiressing  the  data  and  (2)  how  accurate  each  estimator  was. 
Analysis  of  the  residuals  was  done  after  convergence  to  determine  if  the  mean  error  was  approximate^ 
zero  and  if  the  RMS  error  was  close  to  the  one  sigma  values  for  range,  azimuth  and  elevatiotL  Also  of 
inqxtrtance  was  the  accuracy  of  the  estimate  as  measured  by  the  position  and  velocity  errors  relative  to  the 
true  state. 

While  discussing  the  performance  of  the  estimators  it  will  be  useful  to  refer  to  the  first  three 
estimators  ('ets4',  'ets4a',  and  'ets4b')  as  the  fourth  order  estimators.  This  is  because  their  behavim'  tends 
to  be  quite  similar  and  all  three  are  based  on  the  fourth  order  Taylor  Series  iqrptoximation  of  the  two-body 
equations  of  motion.  Likewise  the  last  estimator,  'ets5',  will  sometimes  be  referred  to  as  the  fifth  order 
estimator. 

Because  of  the  non-linear  relationship  between  trade  length  and  maximum  pass  elevation,  it  is 
necessary  to  look  at  the  results  ftom  both  perqrectives.  As  an  example,  trade  lengths  for  the  low  orbit 
used  in  this  study  range  from  6.5-12.4  minutes.  An  estimator  may  give  good  results  for  trades  up  to  1 1.3 
minutes  long  -  which  would  seem  to  cover  a  significant  portion  of  the  possible  passes.  However,  in  terms 
of  pass  elevation  this  would  translate  to  good  performance  only  for  passes  with  a  maximum  elevation  df 
20  degrees  or  less. 

4.2  Estimator  Performance  Using  Perfect  Data. 

Data  rate  proved  to  be  uninqwrtant  when  the  estimators  were  provided  with  noiseless  data.  The  RMS 
and  mean  residual  statistics  along  with  the  position  and  velodty  errors  were  practically  identical 
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regaidless  of  data  rate  and  varied  by  only  a  few  percent.  Because  of  this,  only  information  from  the  low 
data  rate  will  be  presented  and  discussed  for  both  the  low  and  high  orbits. 

4.2. 1  Low  Orbit.  Figures  7  and  8  show  the  RMS  and  mean  error  statistics  of  each  estimator  versus 
track  length  for  the  low  orbit  (satellite  20003).  Figures  9  aitd  10  present  this  same  inftMmation  plotted 
against  the  maximum  elevation  of  the  pass.  As  can  be  seen  from  these  figures,  ’etsS’  consistently 
outperforms  the  other  three  estimators  with  RMS  and  mean  residuals  that  are  approximate^  zero  for  all 
passes.  Although  the  azimuth  and  elevation  residuals  for  the  fourth  order  estimators  diverge  significantly 
d!mcy  fit>m  zero,  they  are  still  relatively  small.  However,  the  range  residuals  indicate  a  definite  weakness 
in  that  parameter.  All  three  of  the  fourth  order  estimators  either  exceed  or  approach  RMS  residuals  of 
100  meters  for  passes  that  last  12.2  minutes  or  more  *-  which  equates  to  elevations  of  48  degrees  or 
greater. 

The  RMS  and  Mean  residuals  indicate  that  the  interval  of  convergence  of  the  fourth  order  estimators 
is  significantly  less  than  most  of  the  track  durations  looked  at  in  the  low  orbit.  Residuals  for  observations 
that  occur  in  a  small  time  interval  (the  interval  of  convergence)  centered  around  the  estimate's  qxich  time 
ate  rather  small  while  residuals  for  observations  outside  this  interval  increase  rapidly  as  At  increases.  As 
the  track  length  increasingly  exceeds  the  interval  of  convergence,  the  percentage  of  the  total  residuals 
outside  the  interval  also  increases  as  does  their  magnitude.  Eventually  these  large  residuals  become 
significant  and  cause  the  RMS  and  Mean  residuals  to  diverge  away  from  zero.  TIk  fifth  order  estimator 
performs  well  because  its  interval  of  convergence  is  significantly  larger.  However,  it  should  be  noted  that 
even  the  fifth  order  estimator’s  residuals  show  a  slight  divergence  away  from  zero.  Even  though  the 
interval  of  convergence  is  larger  it  is  still  finite. 

Figures  1 1  and  12  show  the  position  and  velocity  enors  in  the  estimated  state  from  each  estimator 
versus  track  length  and  maximum  pass  elevation  respectively.  Agaii>,  the  fifth  order  estimator 
significantly  outperforms  the  other  estimators.  Interestingly,  while  the  velocity  error  increases  with  trade 
length  for  'ets4a'  and  'ets4b',  it  decreases  for  'ets4'.  For  tracks  longer  than  1 1  minutes  (about  18  degrees 
in  elevation)  'ets4'  actually  produces  smaller  velocity  errors  than  both  'ets4a'  and  'ets4b'.  The  maximum 
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error  i»oduced  'etsS'  is  1 .8  meters  in  position  and  2.2  cm/sec  in  velocity.  A  more  quantitative  analysis 
of  the  information  contained  in  the  graphs  is  presented  in  Table  9.  Because  cS  the  very  different  shape  ttf 
the  curves  dqwnding  on  whether  the  errors  are  plotted  against  time  or  elevatitm,  weighted  averages  woe 
calculated  for  each  case.  TAle  9  shows  that  the  average  velocity  error  for  'ets4'  is  still  worse  than  any  of 
the  other  estinuitors. 


Table  9.  Position  and  Velocity  Errors  from  Perfect  Data  -  Low  Orbit 


Position  Error  (m 

Velocity  Error  (cm/sec)  I 

Estimator 

Min 

Max 

Avgvs 

Time 

Avgvs 

Elevation 

Min 

Max 

Avgvs 

Time 

Avgvs 

Elevation 

ets4 

46.1 

128.9 

77.8 

80.4 

44.5 

104 

63.5 

53.1 

ets4a 

0.7 

60.2 

18.2 

26.6 

KB 

64.3 

26.4 

33.3 

ets4b 

1.4 

79.9 

24.5 

35.7 

9.8 

89.7 

37.1 

47.2 

etsS 

0.3 

1.8 

0.9 

0^ 

0.3 

2.3 

1.1 

1 

Normally,  it  is  expected  that  as  the  track  duration  increases  the  accuracy  of  the  estimated  state  will 
also  increase  since  more  data  is  available  for  the  estimator.  However,  as  can  be  seen  in  Figure  11,  this  is 
not  true  for  the  fourth  order  estimators.  As  was  the  case  for  the  residuals,  the  interval  of  convergence  for 
these  estimators  is  too  small  and  so  estimates  of  the  state  become  worse  as  the  track  duration  increases. 

The  unexpected  trend  in  velocity  errors  for  'ets4'  relative  to  'ets4a'  and  ’ets4b'  could  possibly  be  related 
to  changes  in  the  value  for  the  J2  potential  function.  Recollect  that  the  J2  potential  function  is  given  by 
Eq  (40).  If  the  potential  at  the  middle  of  the  track  is  considered  indicative  of  the  potential  over  the  entire 
span  of  the  track,  then  a  comparison  of  the  potentials  associated  with  each  pass  will  show  that  at  12 
minutes  (Pass  #7)  the  J2  potential  function  is  a  minimum.  This  could  possibly  reduce  the  inqwrtance  of 
the  J2  terms  contained  in  'ets4a'  and  'ets4b'.  Additionally,  the  J2  terms  could  have  the  effect  of  adding 
numerical  noise  to  the  state  predictions.  The  could  be  a  factor  in  explaining  why  'ets4'  gives  better 
velocity  estimates  than  the  other  fourth  order  estimators  for  passes  longer  than  lO.S  ntinutes. 
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Figure  8.  Converged  Mean  Residuals  vs.  Track  Length  -  Low  Oibit 


Max  Elevation  (deg) 


-B- 

♦ 


ets4 

ets4a 

ets4b 

ets5 


Figure  10.  Converged  Mean  Residuals  vs.  Max  Elevation  --  Low  Oibit 
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Figure  1 1 .  Position  and  Velocity  Errors  vs.  Track  Length  ~  Low  Oibit 
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Figure  12.  Position  and  Velocity  Errors  vs.  Max  Elevation  ~  Low  Orbit 
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4.2.2  High  Orbit.  Figures  13  and  14  s1k>w  the  RMS  and  mean  enor  statistics  each  estimator 
versus  track  length  for  the  high  orbit  (satellite  20(X)4).  Figures  IS  and  16  present  this  same  information 
plotted  against  the  maximum  elevation  of  the  pass.  As  can  be  seen  from  these  figures,  'etsS'  continues  to 
significantly  outperform  the  other  three  estimators. 

The  fourth  order  estimators  show  a  definite  irmbility  to  handle  the  significantly  kmger  trade 
durations.  Elevation  residuals,  both  RMS  and  mean,  .'i^i^roach  the  one  sigma  value  of  0.02S  degrees  at 
elevations  above  60  degrees  (pass  length  of  21  minutes)  and  exceed  0.01  degrees  for  passes  dxive  30 
degrees  (20  minutes).  A  uth  residuals  are  also  very  hi^  but  exhibit  some  interesting  behavior.  The 
RMS  residuals  for  azimuth  have  a  maximum  for  tracks  lasting  about  21  minutes  (max  elevation  of  48 
degrees)  and  an>roach  the  ‘etsS*  value  at  pass  elevations  close  to  90  degrees.  For  passes  between  15-80 
degrees,  'ets4'  actually  has  a  lower  RMS  value  for  azimuth  residuals  than  either  'ets4a'  or  'ets4b'.  The 
mean  residual  values  peak  for  tracks  lasting  about  17.S  minutes  (16  degrees  elevation)  and  pass  through 
zero  for  20  minute  (30  degree)  passes  before  once  again  increasing.  Range  residuals,  both  RMS  and 
mean,  show  a  steady  growth  with  respect  to  time  and  elevation.  The  maximum  RMS  residuals  go  as  high 
as  900  meters  and  the  one  sigma  value  of  100  meters  is  exceeded  for  tracks  lasting  longer  than  17.5 
minutes  (about  16  degrees  elevation). 

Figures  17  and  18  show  the  position  and  velocity  errors  in  the  estimated  state  from  each 
estimator  versus  track  length  and  maximum  pass  elevation  reflectively.  As  expected,  the  fifth  order 
estimator  significantly  outperforms  the  other  estimators.  Again,  'ets4'  had  smaller  velocity  errors  than 
'ets4a'  or  'ets4b'  after  the  track  length  had  increased  to  a  certain  level.  In  this  case,  the  crossover  point 
was  at  a  track  length  of  about  14  minutes  (8  degrees  in  elevation).  As  was  seen  for  the  range  residuals 
earlier,  the  fourth  order  estimators  show  rapid  growth  in  position  and  velocity  errors  to  very  high  levels. 
Position  errors  exceeded  100  meters  for  tracks  of  14  minutes  (about  10  degrees  elevation)  or  more.  The 
information  in  Table  10  provides  a  quantitative  look  at  this  data  for  comparative  purposes.  Because  of  the 
very  different  shape  of  the  curves  depending  on  whether  the  errors  are  plotted  against  time  or  elevation, 
weighted  averages  were  calculated  for  each  case. 
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Table  10.  Position  and  Velocity  Errors  from  Perfect  Data  -  High  Obit 


Position  Error  (m 

Velocity  Error  (cm/sec)  I 

Estimator 

Min 

Max 

Avg  vs 
Time 

Avg  vs 
Elevation 

Min 

Max 

Avg  vs 
Time 

Avg  vs 
Elevation 

ets4 

29.3 

755.4 

251.7 

348.2 

89.1 

380.7 

174.8 

215.9 

ets4a 

wsm 

657.2 

186.4 

285.1 

26.9 

3%.4 

142.3 

190.7 

ets4b 

mm 

722.9 

205.3 

314.2 

29.4 

449.3 

160.6 

215.9 

ets5 

2 

23.6 

9.3 

9.5 

1.5 

21.7 

8.9 

8 

The  same  relationship  between  track  length  and  order  of  approximation  in  the  equations  of  state 
that  was  seen  for  the  low  orbit  is  even  more  evident  in  the  high  orbit  -  both  for  residuals  and  for  accuracy. 
It  should  also  be  noted  that  'etsS'  and  'ets4b'  both  use  the  same  state  transition  matrix  but  'etsS'  uses  a 
higher  order  approximation  for  the  equations  of  state.  Therefore,  the  dramatic  improvement  in 
performance  between  the  two  estimators  demonstrates  that  estimator  performance  is  more  sensitive  to  the 
"accuracy”  of  the  equations  of  state  than  to  the  state  transition  matrix. 
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14.  Converged  Mean  Residuals  vs.  Track  Length  -  High  Orbit 
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Figure  16.  Converged  Mean  Residuals  vs.  Max  Elevcdon  ~  High  Orbit 
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Figure  17.  Position  and  Velocity  Errors  vs.  Track  Length  ~  High  Orbit 
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Velocity  Error 


4.3  Estimator  Performance  Using  Noisy  Data 

In  order  to  ensure  that  the  random  number  generator  had  not  biased  the  data  to  give  eq)eciallv 
good  or  poor  results,  several  data  runs  were  performed  for  each  combination  cX  orbit  and  data  rate  and  the 
results  averaged  together.  Each  data  tun  used  a  different  initial  seed  number  for  the  random  nuttiber 
generator  so  that  a  totally  different  set  trf'  random  numbers  were  generated  during  each  data  run.  Twelve 
data  runs  were  done  for  both  the  high  and  low  orbits  using  the  fast  data  rate.  Five  tuns  were  done  for  each 
orbit  using  the  slow  data  rate. 

Testing  of  the  estimators  using  noise-less  data  indicated  that  only  the  fifth  order  estimator  was 
suitable  for  use.  Because  of  this  the  fourth  order  estimators  were  dropped  from  further  analysis  and  only 
the  fifth  order  estimator  was  analyzed  using  noisy  data. 

It  was  observed  that  when  the  observations  contained  no  noise,  the  data  rate  was  not  a  ftKtor  in 
estimator  performance.  However,  when  the  observations  do  include  noise,  the  higher  data  rate  results  in 
better  estimator  performance. 

4.3. 1  Residuals.  A  typical  plot  of  the  average  mean  and  RMS  residuals  for  range,  azimuth  and 
elevation  versus  maximum  pass  elevation  after  convergence  is  shown  in  Figure  19.  As  can  be  seen,  the 
average  mean  residuals  are  all  very  close  to  zero  and  the  RMS  residuals  are  all  very  close  to  their  expected 
one  sigma  values.  This  behavior,  with  only  minor  differetKes,  was  exhibited  in  all  cases  regardless  of 
orbital  height,  data  rate,  maximum  pass  elevation  or  track  duration. 

Figure  20  shows  a  typical  plot  of  the  RMS  residuals  before  and  after  differential  correction. 
Although,  the  final  pass  RMS  residuals  all  seem  to  be  zero,  they  are  actually  the  same  as  in  Figure  19. 

The  size  of  the  pre-correction  RMS  residuals  dominates  the  vertical  scale  thus  causing  the  post-correction 
values  to  appear  to  be  zero. 

Figures  2 1  and  22  shows  the  actual  pre  and  post  correction  residuals  in  range,  azimuth,  and 
elevation  for  single  passes  of  satellite  20004  over  the  site.  Figure  21  is  a  minimum  elevation  pass 
(elevation  =  3.6  deg)  while  Figure  22  is  a  maximum  elevation  pass  (elevation  =  87.7  deg).  Both  passes 
were  generated  using  the  fast  data  rate.  Again,  these  plots  are  typical  of  all  the  cases  looked  at.  The 
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"spike”  in  the  first  pass  residuals  seen  in  Figure  22  are  caused  by  the  r^d  change  in  azinuith  and 
elevation  as  the  satellite  passes  almost  directly  overhead  of  the  site.  As  can  be  seen,  in  comparison  to  the 
pre-correction  residuals,  the  residuals  after  correction  are  very  small  and  are  a|q>roximateIy  zero. 

Figure  23  plots  the  low  elevation  pass  post-correction  residuals  used  in  Figure  21  versus  their 
e^qiected  one  sigma  limits.  Although  outliers  exist,  the  greater  majority  of  residuals  lie  within  the  one 
signui  limits.  It  should  also  be  noted  that  the  outliers  all  af^rear  to  be  within  three  sigma  of  zero.  This  is 
exactly  what  is  expected  since  after  correction  the  major  source  for  the  residuals  should  be  the  noise  that 
was  in  the  observations. 

All  of  the  information  presented  in  Figures  19-23  indicates  that  the  fifth  order  estimator  has 
effectively  extracted  all  of  the  pertinent  information  from  the  observation  data. 
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Figure  19.  Average  Mean  and  RMS  Residuals  After  Convergence  versus  Max  Elevation 
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Figure  21.  Residual  Values  Before  and  After  Correction  •  Low  El  Pass 


56 


Figure  23.  Residual  Values  versus  One  Sigma  After  Convergence  -  Low  El  Pass 
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4.3.2  Accuracy  of  Estimate.  The  estimator's  accuracy  is  evaluated  by  comparing  the  estimated 
state  for  each  pass  to  the  trtie  state  and  calculating  the  position  and  velocity  error  vector  magnitudes. 
Figure  24  plots  the  position  and  velocity  errors  versus  maximum  elevation  for  the  low  oibit  (satellite 
20(X)3)  for  both  the  high  and  low  data  rates.  As  expected,  the  accuracy  of  the  estimated  state  improves  as 
the  maximum  pass  elevation  increases  and  as  the  data  rate  increases.  The  improvement  is  most  rapid  at 
the  lower  elevations  and  seems  to  level  off  at  the  higher  elevations.  This  is  because  of  the  non-linear 
relationship  between  maximum  elevation  and  track  duration.  Breaking  the  plots  in  Figure  24  into  three 
zones;  Zone  1  accounts  for  an  average  of  58%  of  the  total  change  in  enor  magnitude,  Zone  2  accounts 
for  about  28%  of  the  total  change,  and  Zone  3  accounts  for  the  remaining  14%. 

Figure  25  plots  the  position  and  velocity  errors  versus  maximum  elevation  for  the  high  orbit 
(satellite  20004)  for  both  the  high  and  low  data  rates.  Again,  the  accuracy  of  the  estimate  improves  with 
increasing  elevation  and  data  rate  and  the  improvement  is  still  most  rapid  at  low  elevatiohs. 

When  using  the  low  data  rate,  the  plot  of  error  magnitude  versus  elevation  results  in  fairly 
smooth  curves.  However,  the  higher  data  rate  results  in  curves  with  very  different  behavior.  Accuracy 
improves  rapidly  until  the  elevation  reaches  about  10  degrees.  The  accuracy  then  decreases  until  the 
elevation  reaches  15-20  degrees  after  which  it  again  improves.  This  behavior  is  most  noticeable  in  the 
velocity  errors  and  is  most  pronounced  for  the  high  orbit.  In  general,  the  accuracy  of  the  estimate  is  better 
for  the  low  orbit  than  for  the  high  orbit.  Position  errors  are  consistently  about  10-16  meters  higher  for  the 
high  orbit.  The  velocity  errors  show  a  much  larger  degree  of  variability.  At  elevations  less  than  10 
degrees,  the  velocity  errors  for  the  high  orbit  are  actually  less  than  for  the  low  orbit.  Between  16-35 
degrees,  the  velocity  errors  for  the  high  orbit  are  about  50%  greater  than  for  the  low  orbit.  For  elevations 
above  50  degrees,  the  velocity  errors  for  the  high  and  low  orbits  vary  by  less  than  10%.  Figure  26 
graphically  compares  the  errors  for  the  high  and  low  orbits. 
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Figure  24.  Position  and  Velocity  Errors  versus  Max  Elevation  •  Low  Oibit 
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Figure  25.  Position  and  Velocity  Errors  versus  Max  Elevation  •  High  Orbit 
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Figure  26.  Position  and  Velocity  Errors  versus 


Max  Elevation  •  High  Data  Rate 


4.3.3  Covariance  Matrix.  The  analysis  in  the  previous  section  was  possible  only  because  the 
true  state  was  a  known  quantity  from  the  truth  model.  In  reality,  the  true  state  is  usually  not  known  and 


so  some  other  way  of  assessing  the  quality  of  the  estimate  is  necessary.  This  can  be  found  in  the 
covariance  matrix  which  is  generated  during  the  estimation  process. 

Using  the  state  covariance  matrix  it  is  possible  to  calculate  the  covariance  associated  with  the 
position  and  velocity  errors  (4:43).  The  first  stq>  is  to  write  an  expression  relating  the  error  to  the  slate 
vector.  Next  calculate  the  gradient  of  this  expression  with  respect  to  the  state  vector.  Now  the  gradiem 
and  state  transition  matrix  can  be  combined  to  find  the  variaiKe  in  each  error; 

=V/?’'P^(/o)V/?  (57) 

where  A  is  the  scalar  expression  relating  the  position  or  velocity  error  magnitude  to  the  state  vector  and 
PgtateC^o)  ^  covariance  matrix  associated  with  the  estimate. 

Table  1 1  compares  the  standard  deviation  in  position  and  velocity  (as  predicted  by  the  covariance 
matrix)  with  the  actual  average  errors  for  six  different  pass  profiles.  The  estimator  does  a  good  job  of 
estimating  the  variance  (which  is  merely  the  square  of  the  standard  deviation)  for  position  but  is  very 
qrtimistic  when  estimating  the  variance  for  velocity. 


Table  1 1 .  Comparisoii  of  Predicted  and  Actual  Estimator  Accuracy 
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3.116 
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6.269 

N 

87.78 

3.53 

17.7 

5.014 

0.452 

8.294 

18.350 

The  estimator  is  slightly  optimistic  in  calculating  the  variance  in  position  for  most  passes,  but  as 
the  track  durations  become  very  long  it  seems  to  become  more  and  more  optimistic.  This  can  be  seen  in 
the  last  two  rows  of  Table  11.  This  probably  indicates  that  the  track  lengths  are  beginning  to  exceed  the 
span  over  which  the  approximate  equations  of  state  give  a  good  fit  with  reality. 


Estimates  of  the  variance  for  velocity  are  very  optimistic  and  differ  from  the  actual  errors  by  an 
order  of  magnitude  on  average.  Both  high  and  low  orbits  show  increasingly  optimistic  estirnmes  the 
variance  as  the  track  length  increases.  Interestingly,  the  estimator  generally  does  a  better  job  with  the 
higher  orbit  although  the  effect  of  track  length  mentioned  before  seems  to  be  manifesting  itself  in  the  last 
pass  profile.  The  most  likely  reason  for  the  overly  optimistic  estimate  for  the  velocity  variance  is  that  the 
velocity  terms  in  the  equations  of  state  were  only  developed  to  the  fourth  order  while  the  position  terms 
were  developed  to  the  Mh  order. 

4.4  Summary 

The  analysis  was  begun  with  four  different  estimators  to  test  different  orders  of  approximation  in 
the  equations  of  state  and  state  transition  matrix.  Measures  of  estimator  performance  were  compiled 
using  noise-less  or  "perfect"  observations.  Based  on  these  results,  only  one  of  the  estimators  was  chosen 
for  further  analysis.  When  given  observation  data  that  included  gaussian  noise,  the  fifth  order  estimator 
effectively  extracted  all  of  the  pertinent  information  from  the  data  and  provided  accurate  estimates  of 
position  and  velocity. 
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V.  Conclusiwis  and  Recommendations 


5.1  Introduction 

The  analysis  of  Chapter  4  shows  that  a  diCTeieiUial  corrector  based  on  Taylor  Series  iq)proxiniations  of 
the  equations  of  state  and  state  transition  matrix  can  adequately  estimate  the  state  for  near  earth  orbits. 

The  key  element  in  the  success  of  the  estimator  was  the  use  of  a  sufficient  number  of  terms  in  the 
approximate  equations. 

5.2  Conclusions 

The  number  of  terms  in  the  Taylor  Series  equations  state  required  to  properly  nradel  the  orbital 
dynamics  is  ultimately  driven  by  the  eiqiected  track  duration.  This  is  because  the  Taylor  Series  is  an 
expansion  in  time.  As  the  time  span  over  which  the  equations  must  operate  increases  (i.e.  longer  tracks), 
additional  terms  must  be  added  to  the  equations  or  they  will  fail  to  jHoperly  model  the  dytuunics. 

Although  the  fifth  order  estimator  rqrpears  to  perform  well  for  the  two  orbits  studied,  it  would  prDbid>ly 
begin  to  perform  poorly  as  the  orbital  altitude  increases  since  altitude  is  a  ctor  in  track  length. 

As  expected,  the  two  body  effect  is  the  dominant  term  in  the  state  and  Ui-  x  tranation  matrix 
equations.  However,  although  small,  the  effect  of  the  J2  geqxrtential  term  is  noticeable  for  these  low 
earth  orbits  and  therefore  the  estimator  dynamics  needs  to  take  this  into  account.  This  does  not  mean  that 
the  estimator  will  not  work  without  the  J2  term  in  the  dynamics  -  only  that  it  will  not  work  quite  as  well. 

Estimator  performance  is  much  more  sensitive  to  the  order  of  the  equatiorts  of  state  than  to  the  order 
ci  the  state  transition  matrix.  This  is  demonstrated  by  the  difference  in  performance  between  'etsS'  and 
'ets4b'.  Both  use  the  sariK  state  transition  matrix  equations  but  'etsS'  has  equations  of  state  one  order 
higher  than  those  of  'ets4b'  (see  Table  2  and  section  4.2). 

As  expected,  the  accuracy  of  the  estimator  was  related  to  the  number  of  available  observations.  The 
higher  data  rate  resulted  in  consistently  more  accurate  estimates  than  the  lower  data  rate.  As  the  trade 
length  increased  due  to  increasing  maximum  pass  elevation  the  accuracy  of  the  estimate  also  increased. 
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Accuracy  was  also  affected  I9  orbit  altitude  and  tended  to  be  better  for  the  lower  (»bit.  This  was  {Hobably 
due  to  reduced  effectiveness  of  the  approximate  equations  over  the  longer  time  s^Mns  gmerated  by  the 
higher  orbit. 

Two  effects  can  be  seen  to  be  working  against  each  other.  Longer  trade  lengths  mean  mme  data  is 
available  for  the  estimator  and  so  the  accuracy  of  the  estimate  is  better.  However,  longer  track  lengths 
also  mean  longer  time  spans.  As  noted  before,  the  use  of  approximate  equations  for  the  dynamics  means 
that  the  effectiveness  of  the  dyiuimics  model  decreases  with  increasing  time  span.  So  one  effect  works  to 
inqnove  the  estimate  while  the  other  works  to  degrade  the  estimate. 

As  developed  in  this  study,  the  estimator's  covariance  matrix  provides  overly  optimistic  estimates  of 
the  accuracy  of  the  estimated  velocity.  The  most  likefy  explanation  for  this  is  that  the  state  and  state 
transition  matrix  equations  for  velodty  are  one  ortter  lower  than  the  equations  for  position. 

5.3  Recommendations 

It  should  be  possible  to  improve  the  accuracy  of  the  covariance  matrix  improving  the  velodty 

equations  to  be  of  the  same  order  as  the  position  equations.  This  should  be  done  for  both  the  equations  of 
state  and  the  state  transition  matrix.  Only  additional  terms  for  the  two-body  effea  would  need  to  be  added 
because  of  the  relatively  small  effect  of  J2. 

The  order  of  the  Taylor  Series  approximations  should  probably  be  increased  to  sixth  order  or  possibly 
higher.  This  is  because  the  performance  of  the  fifth  order  estimator  is  expected  to  degrade  when  orbits 
higher  than  those  looked  at  in  this  study  must  be  estimated.  Note  that  the  highest  orbit  in  this  study  has  a 
period  of  1 10  minutes  while  the  definition  of  Near  Earth  includes  satellites  with  periods  up  to  225 
minutes. 

Consider  using  a  numerical  integrator  to  pre^gate  the  state  and  state  transition  matrix  instead  of  a 
Taylor  Series  approximation.  Not  only  would  the  calculation  of  the  state  and  state  transition  matrix  be 
more  exact,  but  the  problem  of  an  effective  interval  of  convergence  associated  with  using  approximate 
equations  would  be  eliminated.  Also,  the  major  reason  for  using  approximate  equations  was  to  reduce 
processing  time  since  numerical  integrators  are  usually  very  computationally  intensive.  However,  as  the 
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order  of  the  approximate  equations  increases,  the  equations  become  quite  com|riex  and  kMe  their  time 
saving  advantage  over  numerical  integrati<m. 

5.4  Potential  Areas  of  Study 

Additional  testing  on  different  orbits  should  be  done.  Higher  altitude  orbits  shoukl  be  looked  at  to 
determine  if  additional  terms  are  needed  in  the  Taylor  Series  approximations  for  the  equations  of  state  and 
state  transition  matrix.  Other  orbital  aspects  such  as  inclination  and  eccentricity  should  also  be 
investigated.  Of  special  interest  would  be  orbits  with  inclinations  of  90  and  63.4  degrees  since  these  orbits 
minimize  the  J2  effect  on  rotation  of  the  line-of-nodes  and  line-of-jq)sides  reqrectively. 

Further  analysis  on  the  effect  of  drag  on  the  estimator's  accuracy  should  also  be  dotte.  This  could  be 
done  by  repeating  the  same  sequence  of  data  runs  done  for  this  study  but  using  truth  data  that  does  not 
itKlude  drag  effects.  The  results  could  then  be  compared  and  the  drag  effects  determined. 

Investigation  of  the  "cup"  seen  in  the  velocity  errors  near  the  10  degree  maximum  pass  elevation 
could  be  done  to  explain  why  this  occurred.  Also,  additional  passes  could  be  analyzed  to  better 
characterize  the  relationship  between  pass  elevation  and  estimator  performance. 

The  estimator  could  be  rewritten  to  use  a  numerical  integrator  instead  of  the  Taylor  Series 
^roximations  for  the  state  and  state  transition  matrix  equations.  Of  primary  interest  would  be  whether 
this  would  improve  the  accuracy  of  the  estimate  and  its  associated  covariance  matrix  as  well  as 
eliminating  the  interval  of  convergence  problem. 

One  of  the  basic  premises  of  this  study  is  that  the  output  from  the  estimator  can  be  groiq)ed  with  other 
similar  estimates  to  solve  for  other  perturbations  such  as  drag.  Develc^ment  of  a  batch  or  sequential 
estimator  to  use  this  data  to  produce  elsets  similar  to  those  currently  produced  by  the  SSC  should  be 
undertaken.  Analysis  of  this  "master"  estimator  should  probably  focus  on;  (1)  obtaining  good  estimates  of 
drag,  (2)  generating  elsets  compatible  with  and  comparable  to  those  produced  by  the  SSC,  (3)  determining 
if  SSN  performance  would  be  improved  by  using  this  approach. 
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Currently  the  SSC  is  not  really  set  up  to  use  a  state  vector  and  covariance  matrix  as  input  to  their 
estimator  (known  as  'SGP4').  Therefore  it  might  be  useful  to  see  if  the  state  vector  could  be  converted  to  a 
form  more  readily  used  by  the  SSC  algorithms. 

5.5  Summary 

This  study  demonstrated  that  a  least  squares  estimator  using  approximate  equations  for  the  state  and 
state  transition  matrix  can  effectively  and  accurately  estimate  the  state  of  a  Near  Earth  satellite  using 
observations  from  a  ground  based  space  surveillance  radar.  Some  limitations  of  the  e^imator  used  in  the 
study  and  ways  to  correct  or  eliminate  them  were  also  identified.  Further  areas  of  study  were  suggested  to 
fill  in  gaps  that  could  not  be  fully  covered  in  this  study.  Of  special  note  was  the  increasing  conq)lexity  in 
the  aiq>roximate  equations  required  to  handle  longer  and  longer  sensor  trades  and  the  potential  of  using  a 
numerical  integrator  to  eliminate  this  problem. 


Appendix  A.  Coordinate  Frame  Tran^ormations 


A.  I  Determining  LMST  and  Site  ECI  Coordinates 

At  any  given  time  the  earth  will  have  a  certain  orientation  relative  to  the  ECI  coordinate  axes.  The 
Greenwich  mean  sidereal  time  (GMST)  qiecifies  the  angle  along  the  celestial  equator  (measured 
eastward)  between  the  mean  equinox  and  the  Greenwich  meridian.  The  GMST  for  any  arbitrary  time  is 

IJT 

0,  =  100.4606184  +  36(X)0.77004J  + 0.000387933 +360.98564724—  (58) 

*  24 


where  6g  is  the  GMST  in  degrees,  J  is  the  number  of  Julian  centuries  from  Julian  year  2000.0,  and  UT  is 
the  universal  time  in  hours  (8:21,  23).  The  Julian  century  can  be  calculated  from 

,  Jo -2451545.0 

J  =  — -  (59) 

36525 

where  Jq  is  the  value  of  the  Julian  date  at  0^  UT  (8:  21).  Once  the  GMST  is  known,  the  site's  local  mean 


sidereal  time  (LMSTi  can  be  calculated  as  6  -  6,.  +  A.  where  A  is  the  site's  east  loneitude.  Fieure  27 
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The  eaith  can  be  modeled  as  an  oblate  q>heroid  with  an  ellipsoidal  cross  section.  The  seini>niajor 
axis  lies  along  the  earth's  equatorial  radius  and  the  semi-minor  axis  is  the  earth's  polar  radius  (2;  93). 
Since  the  cross  section  through  the  equator  is  still  circular  the  longitude  of  a  sensor  site  is  unchanged. 
However,  by  using  this  model  the  determination  of  the  latitude  of  a  site  on  the  earth's  surface  is  mcHe 
complex  and  leads  to  the  definitions  of  geocentric  artd  geodetic  latitudes.  As  can  be  seen  in  Figure  28, 
geocentric  (L')  and  geodetic  (L)  latitudes  are  quite  dififerent.  When  the  latitude  for  a  site  is  given,  it  is 
usually  the  geodetic  latitude  that  is  meant  (2;  94).  The  reference  ellipsoid  used  in  this  model  has  the 
following  dimensions:  Equatorial  radius  (a^)  =  6378. 143  km.  Polar  radius  (b^)  =  6336.785  km,  and 
Eccentricity  (e)  =  0.08182  (2:94). 


If  the  equatorial  axis  is  labeled  "x"  and  the  polar  axis  is  labeled  "z",  then  the  rectangular  coordinates 
of  a  site  will  be 


X  = 


z  = 


sin^  L 
o.(l-e') 


+  h 


yl\-e^  sin^  L 


+  h 


cosL 


sinZ. 


(60) 
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where  is  the  earth's  equatorial  radius,  e  is  the  earth's  eccentricity,  L  is  the  geodetic  latitude,  and  h  is  the 
site's  elevation  dxwe  sea  level  (2:  98).  When  combined  with  the  LMST,  the  site's  location  in  ECI 
coordinates  is 

R  =  X  cosdi  +  X  sin  9J  +  xk  (61) 

A.2  Earth  Centered  Inertial  (ECI)  to  Topocentric  Tranrformation 

Once  the  site's  inertial  coordinates  are  known  the  satellite's  position  relative  to  the  site  is  computed 
as; 

dr  =  (62) 

where  dr  is  the  satellite's  relative  position,  is  the  satellite's  position  in  ECI  coordinates,  and  is 
the  site's  location  in  ECI  coordinates.  The  satellite's  relative  position  is  then  converted  to  the  SEZ  frame 
using; 

p,  =  xsinXcos0+>'sinXsin0-xcosA, 

p,  =>'cos0-xsin0  (63) 

p,  =  xcosX,cos0+^cosXsin0  +  xsinX 

where  Pg,  die  components  of  the  relative  position  vector  in  the  SEZ  frame,  x,  z  are  the 

components  of  the  relative  position  vector  found  in  Eq  (61),  X.  is  the  site's  latitude,  and  6  is  the  LMST  for 
the  site  at  the  time  of  the  observation.  The  SEZ  vector  is  finally  converted  to  the  topocentric  coordinates 
by; 

range  =  p  =  -^x^  +y^  +z^ 

elevation  =  3  =  asin(p,  /  p)  (64) 

azimuth  =  aco 
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A.3  Topocentric  to  ECl  Transformation 

Convening  from  u^xxxnthc  coordinates  to  ECi  coordinates  is  basically  the  reverse  of  going  from  ECI 
to  topocentric.  The  first  step  is  still  to  find  the  site's  inertial  location.  Once  that  is  dcme  the  topocentric 
coordinates  are  transformed  into  the  SEZ  frame: 


p,  =  -pcosPcosa 

p,  =pcos3sina  (65) 

p,  =  r>sin0 

The  nest  step  is  to  convert  from  SEZ  coordinates  to  ECI  coordina*.  s  using: 

dx  =  p^  sin  A.cos6  -  p,  sin  6  +  p,  cosXcosG 

dy  =  sin  X  sin  0  +  p,  cos0  +  p^  cosX  sin  0  (66) 

<fe  =  p,  cosX  +  p,  sinX. 

where  dr,  dy,  and  d?  are  the  components  for  the  relative  position  vector  between  the  site  and  the  satellite 
expressed  in  ECI  coordinates.  The  final  step  then  is  to  add  the  relative  position  vector  and  the  site 
location  together: 

(67) 
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Appendix  B.  The  Data  Linearization  Matrix,  H 


The  data  linearization  matrix  is  given  by  Eq  (18)  as: 

For  a  satellite  orbiting  the  earth,  the  observation  relation  G  is; 


(18) 


Vp'+P'+Pr 

G  = 

02 

= 

Ji-atan(p,/pj 

[eij 

atan(p7Vp'  +p;)^ 

where  p,  az,  el  are  the  range,  azimuth  and  elevation  components  of  the  t<^xx:entric  observation  vector  and 
Ps>  Pe>  Pz  ^  ^  components  of  the  topocentric  observation  veaor.  In  the  SEZ  topocentric  frame, 
the  satellites  position  is; 


PsEZ 

where  r  is  the  satellite's  position  in  ECI  coordinates,  R  is  the  site's  position  in  ECI  coordinates  and  D  is 
the  rotation  matrix  from  the  ECI  ftame  to  the  SEZ  frame; 


D  = 


sinXcos6 

-sin0 

cosX.cos6 


sin  X  sin  6 

COS0 

cosXsin0 


-cosX 

0 

sinX 


where  A.  is  the  site's  latitude  and  6  is  the  site's  LMST. 


(70) 


Since  G  is  linearized  with  respect  to  the  state  vector,  the  site  position  can  be  ignored  and  so 

6p^  =  DBr 


(71) 


and 


dz  = 


^PsEZ 


=  H,D8r 


(72) 
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where  z  is  the  observation  vector  range,  azimuth  and  elevation.  Changes  in  the  observation  vector  are 

linearly  related  to  changes  in  the  satellites  state  as: 

6i  =  B5»  (73) 

and  so  Eq  (72)  can  be  rewritten  and  solved  for  H  as; 

H  =  H,D  (74) 

since  6r/5x  =  1  . 

After  differentiating  the  equations  in  Eq  (68)  the  result  is: 
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Appendix  C.  State  Transition  Matrix,  <I> 


The  state  transition  matrix,  commonly  denoted  by  <I>,  is  found  by  taking  the  gradient  of  the  equations 
of  state  evaluated  at  the  initial  conditions.  When  the  equations  of  state  are  available  in  e)q)licit  form  this 
relationship  can  be  expressed  as. 

=  (14) 

where  xCtg)  is  the  state  at  the  initial  or  reference  conditions  and  h  is  the  relationship  that  relates  the  state 
at  Iq  to  the  state  at  any  other  time  t. 

In  this  study  the  equations  of  state  are  explicit  expressions  obtained  from  a  Taylor  Series  expansion  of 
the  equations  of  motion.  As  can  be  seen  in  Section  3.3.2,  Eqs  (48,  S1-S6),  these  equations  are  rather  large 
and  complex,  so  the  easiest  way  to  derive  the  <I>  matrix  is  to  develq)  it  term  term  and  then  add  all  the 
parts  together.  For  notational  purposes,  each  term  in  Uk  development  will  be  denoted  as  ,  where  the 
subscript  refers  to  the  order  of  the  term.  The  same  notation  is  used  when  referencing  a  particular  term  in 
the  Taylor  Series  (such  as  bo  which  would  be  the  zero  term  in  the  series).  Also,  all  quantities  are 
evaluated  at  the  reference  or  initial  conditions. 

C.  /  The  Two-body  Terms 
The  zero  order  term  is; 


<Do=V.ho=V|  1=1 


(76) 


where  r,  v  are  the  initial  position  and  velocity  vectors  and  I  is  a  6x6  identity  matrix. 
The  first  order  term  is; 


(77) 
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where  a  is  the  initial  acceleration  vector  (reference  Eq  (53)),  0  is  a  3x3  null  matrix,  I  is  a  3x3  identity 
matrix  and  A  is  the  3x3  matrix; 


(Sx^/r’ -  l/r^)  3xy/r*  3xz/r’ 

A  =  |i  3xylr^  (3y/r*-l/r^)  3^/r*  (78) 

3xr/r*  3>«/r*  (3z7r*-l/r*) 

where  ^  is  the  earth's  gravitational  parameter,  x  y.  z  are  the  position  coordinates  of  the  state  in  the  ECI 
frame  of  reference,  and  r  is  the  magnitude  of  the  position  vector. 

The  second  order  term  is: 


where  a-dot  is  the  first  time  derivative  of  the  acceleration  vector  (see  Eq  (SI)),  the  3x3  matrix  B  is; 


'[l5(rv)jfV»'’-3(rv)-6i«]/r’  [l5(r- v)jiv/r’ -3var  -  3iiy]/r’  [l5(rv)«/r’ -3ti'jr-3ia]/r’ 

B  =  u  [l5(r- v)j9'/r’ -3»w -3ifv]/r’  [l5(r- v)>V'’*  "3(r  v)-6»5']/r‘  [l5(r- v)>«/r* -3ny -3w]/r’ 

[l5(r- v)jB/r’  -3»».*  -3ii*]/r’  [l5(r- v)>x/r’  -3iiy  -3w]/r’  [l5(r •  v)rV'’’  -3(r  v)-6i«]/r’ 

(80) 

and  the  3x3  matrix  C  is: 


-(l/r^  -  3xV»"0  3xy/r*  3xr/r* 

C  =  n  3xy/r^  -3y^  jr^)  ^yzjr^  (81) 

3xr/r*  3W'’*  -(l/r^  -  3zV'’0 

The  third  order  term  is; 


where  ii  is  given  by  Eq  (52)  and  the  3x3  matrix  D  is; 
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D(l,l)=  -\Q5{rdv)^  ! r*  -3jur*/r* +[l5(r<ft')*  +60(rdv)ux  +  l5x^  (veJ^  -ji/r)]/r’ 

-n/r*  -[6«’  +3(ve/^  -H/r)]/r’  jr* 

D(l,2)=  [30(rdv)(vx  +  uy)  +  lSxy(vel^  -H/r)yr’  -6uvlr^  -9\ixylr*  -  105(«A')* ly 
D(l,3)  =  [30(r</w)(wx  +«2)  +  15xz(ve/*  -(i/r)jyr’  -tuwjr^  -9\tx2fr*  -lOSirdv)^  xz 
D(2.1)=  D(l.2) 

D(2,2)=  -105(r<A')’>'’/r’  -3n>'’/r* +[l5(r<A')’  +60(r<A')yy +  15>'’ (ve/’ 

-n/r‘  -[6v^  +3(ve/^  -|A/r)]/r’  +6jiy /r* 

D(2,3)=  [30(rrfv')(Hy +  vz)  +  15>z(ve/^  -^/r)yr’  -6vwfr*  -9n>??/r*  - l05(r<A>)’ yz 
D(3,l)=  D(l,3) 

D(3,2)=  D(2,3) 

D(3,3)=  -lOS(rdv)^  2^  f -3|iz^/r* +[l5(rift')’  +60(ri/w)wz +  15z’ (ve/^  -|i/r)yr’ 

-n/r‘  -[6w’  +3(ve/^  +6\i2^ /r* 

(83) 


where  vel  =  |v|  and  rdv  =  (r  ■  v) . 
The  E  matrix  is  given  by: 


E  =  - 


30(rv)x7'-*-6(rv)-12ttx 
30(r  v)jy/r*  -6vx-6«>' 
30(r- v)xz/r’  -6wx-6uz 


30(r  v)xy/r*  -6vx-6«>' 
30(rv)y7r’-6(rv)-12yy 
30(r-v)yz/r*  -6My-6vr 


30(r  •v)xz/r*  -Swx-eKZ 
30(r  v)yr/r*  -6Hy-6vr 
30(rv)x7r*-6(rv)-12w 


(84) 


The  fourth  order  term  (which  was  only  calculated  for  position)  is: 


C.2  The  J 2  Terms 

Following  the  same  procedure  used  for  the  two-body  elements  of  the  O  matrix,  the  first  order  term  for 
J2  effects  is: 


<1> 


(86) 
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where  >20  ^  ^  ^  ^  the  3x3  matrix  FI  is: 


3/r’  -  15(x’  +z’  )/r’  +  lOSx’z*/''* 

-XSxyfr'  +\(i$xyz^  fr* 

^5xr/r’ +105xz7r* 

Fl  =  -/t 

-15;O»/r’+105j9«Vr* 

3/r’  - 15(>*  +z’  )/r’  +  lOSy'z’/r’ 

-A^ytjr^  +105j«’/r* 

-45xz/r’  +105xzVr’ 

-A^yxfr'  +105>z’/r* 

9lr'-Viz^lr^  +105z7r’ 

(«7) 


The  second  order  term  is; 


(D 


12 


m 


where  ^  is  the  quantity  defined  in  Eq  (39),  *20*^  ^  ^  the  3x3  matrix  F2  is; 


F2(l,l)  =  [-30ux-2w2 -irdv)]/r' -9ASx2^(rdv)/r" 
+105[«r’  +2wx^z+ux2^  +x^  (r<A')]/r’ 
F2(l,2)=  -15(vx  +  «v)/r^  -945><z*(rA-)/r'* 

+  105[v2*  +y{2wxz  +«?’  ■¥x(rdv))yr^ 
F2(l,3)=  -45(i«  +  wx)/r’  -945x"(r<*')/r" 

105r[H7*  +3z(rtft')  +  2wxz  +«r^  +x(r£A’)]/r’ 
F2(2,l)=  -ISCvx  +  M^)^’  -945xr’(r<*')/r" 

+  105^«z’  +y{2wyz  +vz^  +y(rdv)^r^ 
F2(2,2)=  [-30iy  - 2 wz -(«#>')]//’’  -9A5yz^ (rdv)lr" 
+I05[vr^+2»f>’z+vj«^  +y  (r(A')]/r’ 
F2(2.3)  =  -A5(vz  +  wy)lr''  -945z"  {rdv)lr" 

105z[>w^  +3z(rdv)  +  2wyz  +vz*  +>'(r<A')]/r* 
F2(3, 1)  =  -45(«z  +  wx )/r'  - 945xz’  (r<A>)/r" 

105[«z’  +  3x(  wz  ’  +z(rdv))yr* 

F2(3,2)=  -45(vz  +  H’>')/r’  -945,vz’(r</»')/r'‘ 
lOSjvz’  +3>'(m!z*  +  z(r(A'))]/r’ 

F2(3,3)=  -[l80H7+45(r<A')]/r’  -945z*  (r<A')/r" 
210[2i«’+3z’(r<A')]/r’ 
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using  the  same  notation  as  Eq  (83).  The  3x3  F3  matrix  is: 


F3  =  -/l 


3/r’  -  15(*’  +z*)/r’  +105joV 

-15jo»/r’  +105«7'‘* 

-45«/r’+105«Vr’ 


+\05yz^  /r* 

3/r’-15(y’+z')/r’+105>«7r* 
-ASyzy  +  105>«7r* 


-ASxzjr''  +  I05z7»-’ 
-ASytfr''  +  105z7»-’ 
9/r’  -  90z7'-"  +  105zV»-* 


(90) 


The  third  order  term  is; 


The  final  <I>  matrix  is  then; 

<D(rfO  =  Oq  +  <D,  +  Oj  +  Oj  +  O4  +  <I>„  +  <D„  +  Ojj  (92) 

As  was  done  for  the  equations  of  state,  the  equations  for  position  are  essentially  one  order  higher  than 
the  velocity  equations.  Of  special  interest  is  the  fact  that  as  higher  order  terms  are  added  to  the  matrix 
their  derivatives  become  very  large  and  complex  and  increasingly  difficult  to  derive. 
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